# -*- coding: utf-8 -*-
"""
Created on Fri Jan 24 12:02:58 2020

@author: MPecchi
"""

import numpy as np  # make sure you use numpy=1.22 to avoid conflict with shapely=1.7
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pathlib as plib
import itertools as itr
import geopandas as gpd
from scipy.optimize import curve_fit

# =============================================================================
# # function definition: all funciton required in the script are defined here
# =============================================================================
def PathsCreate(file, subfolder=''):
    ''' This function creates 2 folder paths. For each path checks the folder
    existence, if False it creates the folder.

    file : str
        the command __file__ in the script provides the path to the script
        itself. This function builds path starting from the script folder.
    subfolder : add one subfolder level

    in_path : pathlib object
        path to the Input folder.
    out_path : pathlib object
        path to the Output folder.'''
    # start from script folder
    script_path = plib.Path(file).parents[0]  # script folder path
    # create all necessary paths to subfolders starting from the script one
    in_path = plib.Path(script_path, "Input", subfolder)
    out_path = plib.Path(script_path, 'Output', subfolder)
    # check existence of each folders and create them if missing
    plib.Path(in_path).mkdir(parents=True, exist_ok=True)
    plib.Path(out_path).mkdir(parents=True, exist_ok=True)
    return in_path, out_path  # returns the two object-paths


def FigCreate(rows=1, cols=1, plot_type=0, paper_col=1,
              gridspec_hgt_fr=None,
              gridspec_wdt_fr=None, hgt_mltp=1, font='DejaVu Sans'):
    """
    This function creates all the necessary objects to produce plots with
    replicable characteristics.

    Parameters
    ----------
    rows : int, optional
        number of plot rows in the grid. The default is 1.
    cols : int, optional
        number of plot colums in the grid. The default is 1.
    plot_type : int, optional
        one of the different plots available. The default is 0.
        Plot types and their labels:
        0. Std: standard plot (single or grid rows x cols)
        1. Twin-x: secondary axis plot (single or grid rows x cols)
        2. Cmap: colormap plot (single or grid rows x cols)
        3. Ir: Infra-red curves plot (no brokenXaxis)
        4. Ir: Infra-red curves plot (WITH brokenXaxis)
    paper_col : int, optional
        single or double column size for the plot, meaning the actual space
        it will fit in a paper. The default is 1.
    gridspec_wdt_ratio : list of float, optional
        for multiple cols, list of relative width of each subplot.
        The default is None.
    no_par : bool, optional
        if True, no size setting parameters are passed. The default is None.
    font: str
        if the str'Times' is given, it sets Times New Roman as the default
        font for the plot, otherwise the default Dejavu Sans is maintained.
        Default is 'Dejavu Sans'
    hgt_mltp: float
        multiplies the fig height. default is to 1. best using values between
        .65 and 2. may not work with multiplot and paper_col=1 or out of the
        specified range
    Returns
    -------
    fig : object
        the figure object to be passed to FigSave.
    lst_ax : list of axis
        list of axis (it is a list even with 1 axis) on which to plot.
    lst_axt : list of axis
        list of secondary axis (it is a list even with 1 axis)
    fig_par : list of float
        ist of parameters to reserve space around the plot canvas.

    """

    if font == 'Times':  # set Times New Roman as the plot font fot text
        sns.set("paper", font_scale=1.25)
        sns.set_palette('deep', 10)
        sns.set_style("ticks", {'font.family': 'Times New Roman'})
    else:  # leave Dejavu Sans (default) as the plot font fot text
        sns.set("paper", font_scale=1.25)
        sns.set_palette('deep', 10)
        sns.set_style("ticks")
    # single or double column in paperthat the figure will occupy
    if cols > 2:  # numer of columns (thus of plots in the figure)
        raise ValueError('\n FigCreate: cols>2 not supported')

    # width of the figure in inches, it's fixed to keep the same text size
    # is 6, 9, 12 for 1, 1.5, and 3 paper_col (columns in paper)
    fig_wdt = 6*paper_col  # width of the plot in inches
    fig_hgt = 4*paper_col*rows/cols*hgt_mltp  # heigth of the figure in inches
    px = 0.06*(6/fig_wdt)*cols  # set px so that (A) fits the square
    py = px*fig_wdt/fig_hgt/cols*rows/hgt_mltp  # set py so that (A) fits
    # if more rows are added, it increases, but if cols areadded it decreases
    # to maintain the plot ratio
    # set plot margins
    sp_lab_wdt = 0.156/paper_col  # hor. space for labels
    sp_nar_wdt = 0.02294/paper_col  # space narrow no labels (horiz)
    sp_lab_hgt = 0.147/paper_col/rows*cols/hgt_mltp  # space for labels (vert)
    sp_nar_hgt = 0.02/paper_col/rows*cols/hgt_mltp  # space narrow no labels
    # (vert)
    # =========================================================================
    # # 0. Std: standard plot (single or grid rows x cols)
    # =========================================================================
    if plot_type == 0:
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = (0.26*paper_col**2 - 1.09*paper_col + 1.35)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .2/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    # =========================================================================
    # # 1. Twin-x: secondary axis plot (single or grid rows x cols)
    # =========================================================================
    elif plot_type == 1:
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
            lst_axt = [ax.twinx()]  # create a list with secondary axis object
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
            # create list of secondary twin axis
            lst_axt = [axs.twinx() for axs in ax.flatten()]
        # horizontal space between plot in percentage !!! needs DEBUG
        sp_btp_wdt = 1.36*paper_col**2 - 5.28*paper_col + 5.57
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .2/paper_col*cols/hgt_mltp
        # left, bottom, right(DIFFERENT FROM STD), top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_lab_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    # =========================================================================
    # # 2. Cmap: colormap plot
    # =========================================================================
    elif plot_type == 2:
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
            fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_lab_wdt/2, 1-sp_nar_hgt,
                       0, 0, px, py]
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
            fig_par = [None, None, None, None, None, None, px, py]
        lst_axt = None  # no secondary axis in this plot_type
    # =========================================================================
    # # 3. Ir: Infra-red curves plot (no brokenXaxis)
    # =========================================================================
    elif plot_type == 3:
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            raise ValueError('\n FigCreate: plot_type=3: '
                             + 'multiplot not supported')
        lst_axt = None  # no secondary axis in this plot_type
        # left(DIFFERENT FROM STD), bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt/3, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   0, 0, px, py]
    # =========================================================================
    # # 4. Ir: Infra-red curves plot (WITH brokenXaxis)
    # =========================================================================
    elif plot_type == 4:
        # creates two subplots, plots a part of the x-axis in each of them
        # then removes the central splines of the two plots to make it one
        fig, ax = plt.subplots(1, 2, figsize=(fig_wdt, fig_hgt),
                               gridspec_kw={'width_ratios':
                                            gridspec_wdt_fr}, sharey=True)
        lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # set plot margins
        # horizontal space between plot in percentage
        sp_btp_wdt = .05
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt/3, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, 0, px, py]
    elif plot_type == 5:  # plots with different heights
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt),
                               gridspec_kw={'height_ratios': gridspec_hgt_fr})
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = (0.26*paper_col**2 - 1.09*paper_col + 1.35)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .2/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    elif plot_type == 6:  # multiplot without internal x and y tickslabels
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = (0.12*paper_col**2 - 1.09*paper_col + 1.35)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .065/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    elif plot_type == 7:  # multiplot without internal x tickslabels
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt))
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = (0.24*paper_col**2 - 1.09*paper_col + 1.35)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .05/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    elif plot_type == 8:  # plots with different heights
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt),
                               gridspec_kw={'height_ratios': gridspec_hgt_fr})
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = (0.25*paper_col**2 - 1.09*paper_col + 1.35)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .05/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    elif plot_type == 9:  # test for colorbars
        fig, ax = plt.subplots(rows, cols, figsize=(fig_wdt, fig_hgt),
                               gridspec_kw={'width_ratios': gridspec_wdt_fr})
        if rows*cols == 1:  # only 1 plot
            lst_ax = [ax]  # create ax list for uniform iterations over 1 obj.
        elif rows*cols > 1:  # more than one plot
            lst_ax = [axs for axs in ax.flatten()]  # create list of axis
        lst_axt = None  # no secondary axis in this plot_type
        # horizontal space between plot in percentage
        sp_btp_wdt = 0  # (1.26*paper_col**2 - .09*paper_col + 1.55)
        # vertical space between plot in percentage !!! needs DEBUG
        sp_btp_hgt = .3/paper_col*cols/hgt_mltp
        # left, bottom, right, top, widthspace, heightspace
        fig_par = [sp_lab_wdt, sp_lab_hgt, 1-sp_nar_wdt, 1-sp_nar_hgt,
                   sp_btp_wdt, sp_btp_hgt, px, py]
    return fig, lst_ax, lst_axt, fig_par


def FigSave(fig_name, out_path, fig, lst_ax, lst_axt, fig_par,
            xLab=None, yLab=None, ytLab=None,
            xLim=None, yLim=None, ytLim=None,
            xTicks=None, yTicks=None, ytTicks=None,
            xTickLabels=None, yTickLabels=None, ytTickLabels=None,
            legend=None, ncol_leg=1, 
            annotate_lttrs=False, annotate_lttrs_loc='down',
            pdf=False, svg=False, transparency=False,
            subfolder=None, tight_layout=False, title=False,
            ):
    '''
    FIXES:
        1. px, py moved to FIgCreate

    This function takes the obects created in FigCreate and allows to modify
    their appeareance and saving the results.

    Parameters
    ----------
    fig_name : str
        name of figure. It is the name of the png od pfd file to be saved
    out_path : pathlib.Path object. path to the output folder.
    fig : figure object. created in FigSave.
    lst_ax : list of axis. Created in FigCreate
    lst_axt : list of twin (secondary) axis. Created in FigCreate
    fig_par : list of figure parameters for space settings
        left, bottom, right, top, widthspace, heightspace, px, py.
        Created in FigCreate
    tight_layout : bool
        If True, ignore fig_par[0:6] and fit the figure to the tightest layout
        possible. Avoids to lose part of figure, but loses control of margins
    xLab : str.list, optional
        label of the x axis. The default is None.
        can be given as
        0. xLab=None: no axis gets an xlabel
        1. xLab='label': only one str, all axis get the same xlabel
        2. xLab=['label1', None, Label2, ...]: the list must have the size of
            lst_ax and contain labels and or None values. Each axis is
            assigned its label, where None is given, no label is set.
    yLab : str, optional
        label of the y axis. The default is None. Same options as xLab
    ytLab : str, optional
        label of the secondary y-axis. The default is None.
        Same options as xLab
    xLim : list of two values, list of lists, optional
        limits of x axis. The default is None.
        can be given as
        0. xLim=None: no axis gets a xlim
        1. xLab=[a,b]: all axis get the same xlim
        2. xLab=[[a,b], None, [c,d], ...]: the list must have the size of
            lst_ax and contain [a,b] and or None values. Each axis is
            assigned its limit, where None is given, no llimit is set.
    yLim : list of two values, optional
        limits of y axis. The default is None. Same options as xLim
    ytLim : list of two values, optional
        limits of secondary y axis. The default is None.
        Same options as xLim
    xTicks : list of int or float, optional
        list of tiks value to be shown on the axis. The default is None.
    yTicks : list of int or float, optional
        list of tiks value to be shown on the axis. The default is None.
    ytTicks : TYPE, optional
        list of tiks value to be shown on the axis. The default is None.
    legend : str, optional
        contains info on the legend location. To avoid printing the legend
        (also in case it is empty) set it to None.
        The default is 'best'.
    ncol_leg : int, optional
        number of columns in the legend. The default is 1.
    annotate_lttrs : bool, optional
        if True, each plot is assigned a letter between () in the lower left
        corner. The default is False. If a string is given, the string is used
        as the letter in the plot even for single plots.
    annotate_lttrs_loc: str.
        default is 'down', if 'up' is given, the letters are placed on the left
        top corner.
    pdf : bool, optional
        if True, the figure is saved also in pdf in the output folder.
        The default is False, so only a png file with 300dpi is saved
    transparency : bool, optional
        if True, background of PNG figure is transparent, defautls is False.
    subfolder : str, optional
        name of the subfolder inside the output folder where the output will
        be saved. If the folder does not exists, it is created.
        The default is None.
    '''

    fig_adj_par = fig_par[0:6]
    if not any(fig_par[0:6]):  # True if all element in fig_par[0:6] are False
        tight_layout = True
    px = fig_par[6]
    py = fig_par[7]
    n_ax = len(lst_ax)  # number of ax objects
    # for xLab, yLab, ytLab creates a list with same length as n_ax.
    # only one value is given all axis are given the same label
    # if a list is given, each axis is given a different value, where False
    # is specified, no value is given to that particular axis
    vrbls = [xLab, yLab, ytLab, legend]  # collect variables for iteration
    lst_xLab, lst_yLab, lst_ytLab, lst_legend \
        = [], [], [], []  # create lists for iteration
    lst_vrbls = [lst_xLab, lst_yLab, lst_ytLab, lst_legend]  # collect lists
    for vrbl, lst_vrbl in zip(vrbls, lst_vrbls):
        if vrbl is None:  # label is not given for any axis
            lst_vrbl[:] = [None]*n_ax
        else:  # label is given
            if np.size(vrbl) == 1:  # only one value is given
                if type(vrbl) == str:  # create a list before replicating it
                    lst_vrbl[:] = [vrbl]*n_ax  # each axis gets same label
                elif type(vrbl) == list:  # replicate the list
                    lst_vrbl[:] = vrbl*n_ax  # each axis gets same label
            elif np.size(vrbl) == n_ax:  # each axis has been assigned its lab
                lst_vrbl[:] = vrbl  # copy the label inside the list
            else:
                print(vrbl)
                print('Labels/legend size does not match axes number')
    # for xLim, yLim, ytLim creates a list with same length as n_ax.
    # If one list like [a,b] is given, all axis have the same limits, if a list
    # of the same length of the axis is given, each axis has its lim. Where
    # None is given, no lim is set on that axis
    vrbls = [xLim, yLim, ytLim, xTicks, yTicks, ytTicks, xTickLabels,
             yTickLabels, ytTickLabels]  # collect variables for iteration
    lst_xLim, lst_yLim, lst_ytLim, lst_xTicks, lst_yTicks, lst_ytTicks, \
        lst_xTickLabels, lst_yTickLabels, lst_ytTickLabels = \
            [], [], [], [], [], [], [], [], [] # create lists for iteration
    lst_vrbls = [lst_xLim, lst_yLim, lst_ytLim, lst_xTicks, lst_yTicks,
                 lst_ytTicks, lst_xTickLabels, lst_yTickLabels,
                 lst_ytTickLabels]  # collect lists
    for vrbl, lst_vrbl in zip(vrbls, lst_vrbls):
        if vrbl is None:  # limit is not given for any axis
            lst_vrbl[:] = [None]*n_ax
        else:
            # if only list and None are in vrbl, it is [[], None, [], ..]
            # each axis has been assigned its limits
            if any([isinstance(v, (int, float, np.int32, str))
                    for v in vrbl]):
                temporary = []  # necessary to allow append on [:]
                for i in range(n_ax):
                    temporary.append(vrbl)  # give it to all axis
                lst_vrbl[:] = temporary
            else:  # xLim=[[a,b], None, ...] = [list, bool] # no float
                lst_vrbl[:] = vrbl  # a lim for each axis is already given
                
    
    # loops over each axs in the ax array and set the different properties
    for i, axs in enumerate(lst_ax):
        # for each property, if the variable is not false, it is set
        if lst_xLab[i] is not None:
            axs.set_xlabel(lst_xLab[i])
        if lst_yLab[i] is not None:
            axs.set_ylabel(lst_yLab[i])
        if lst_xLim[i] is not None:
            axs.set_xlim([lst_xLim[i][0]*(1 + px) - px*lst_xLim[i][1],
                          lst_xLim[i][1]*(1 + px) - px*lst_xLim[i][0]])
        if lst_yLim[i] is not None:
            axs.set_ylim([lst_yLim[i][0]*(1 + py) - py*lst_yLim[i][1],
                          lst_yLim[i][1]*(1 + py) - py*lst_yLim[i][0]])
        if lst_xTicks[i] is not None:
            axs.set_xticks(lst_xTicks[i])
        if lst_yTicks[i] is not None:
            axs.set_yticks(lst_yTicks[i])
        if lst_xTickLabels[i] is not None:
            axs.set_xticklabels(lst_xTickLabels[i])
        if lst_yTickLabels[i] is not None:
            axs.set_yticklabels(lst_yTickLabels[i])
        if annotate_lttrs is not False:
            if annotate_lttrs_loc == 'down':
                y_lttrs = (py/px*.02)
            elif annotate_lttrs_loc == 'up':
                y_lttrs = 1 - py
            if n_ax == 1:  # if only one plot is given, do not put the letters
                axs.annotate('(' + annotate_lttrs + ')',
                              xycoords='axes fraction',
                              xy=(0, 0), rotation=0, size='large',
                              xytext=(0, y_lttrs), weight='bold')
            elif n_ax > 1:  # if only one plot is given, do not put the letters
                axs.annotate('(' + lttrs[i] + ')', xycoords='axes fraction',
                              xy=(0, 0), rotation=0, size='large',
                              xytext=(0, y_lttrs), weight='bold')
    # if secondary (twin) axis are given, set thier properties
    if lst_axt is not None:
        for i, axst in enumerate(lst_axt):
            axst.grid(False)  # grid is always false on secondaty axis
            # for each property, if the variable is not false, it is set
            if lst_ytLab[i] is not None:
                axst.set_ylabel(lst_ytLab[i])
            if lst_ytLim[i] is not None:
                axst.set_ylim([lst_ytLim[i][0]*(1 + py) - py*lst_ytLim[i][1],
                              lst_ytLim[i][1]*(1 + py) - py*lst_ytLim[i][0]])
            if lst_ytTicks[i] is not None:
                axst.set_yticks(lst_ytTicks[i])
            if lst_ytTickLabels[i] is not None:
                axst.set_yticklabels(lst_ytTickLabels[i])
    # create a legend merging the entries for each couple of ax and axt
    if any(lst_legend):
        if lst_axt is None:  # with no axt, only axs in ax needs a legend
            for i, axs in enumerate(lst_ax):
                axs.legend(loc=lst_legend[i], ncol=ncol_leg)
        else:  # merge the legend for each couple of ax and axt
            i = 0
            for axs, axst in zip(lst_ax, lst_axt):
                hnd_ax, lab_ax = axs.get_legend_handles_labels()
                hnd_axt, lab_axt = axst.get_legend_handles_labels()
                axs.legend(hnd_ax + hnd_axt, lab_ax + lab_axt, loc=lst_legend[i],
                           ncol=ncol_leg)
                i += 1
    try:
        fig.align_labels()  # align labels of subplots, needed only for multi plot
    except AttributeError:
        print('align_labels not performed')
    # if a subfolder is specified, create the subfolder inside the output
    # folder if not already there and save the figure in it
    if subfolder is not None:
        out_path = plib.Path(out_path, subfolder)  # update out_path
        plib.Path(out_path).mkdir(parents=True, exist_ok=True)  # check if
        # folder is there, if not create it
    # set figure margins and save the figure in the output folder
    if tight_layout is False:  # if margins are given sets margins and save
        fig.subplots_adjust(*fig_adj_par[0:6])  # set margins
        plt.savefig(plib.Path(out_path, fig_name + '.png'), dpi=300,
                    transparent=transparency)
        if pdf is not False:  # save also as pdf
            plt.savefig(plib.Path(out_path, fig_name + '.pdf'))
        if svg is not False:  # save also as pdf
            plt.savefig(plib.Path(out_path, fig_name + '.svg'))
    else:  # margins are not given, use a tight layout option and save
        plt.savefig(plib.Path(out_path, fig_name + '.png'),
                    bbox_inches="tight", dpi=300, transparent=transparency)
        if pdf is not False:  # save also as pdf
            plt.savefig(plib.Path(out_path, fig_name + '.pdf'),
                        bbox_inches="tight")
        if svg is not False:  # save also as pdf
            plt.savefig(plib.Path(out_path, fig_name + '.svg'),
                        bbox_inches="tight")
    # add the title after saving, so it's only visible in the console
    if title is True:
        lst_ax[0].annotate(fig_name, xycoords='axes fraction', size='small',
                            xy=(0, 0), xytext=(0.05, .95), clip_on=True)
    return


def latlongdist(lat1, long1, lat2, long2):
    # calculates the distance [km] between two points,
    # use Haversine formula to compute distance between two coordinates
    # (https://www.movable-type.co.uk/scripts/latlong.html)
    R_earth = 6371  # earth radius [km]
    lat1_rad = lat1*np.pi/180  # deg to rad
    long1_rad = long1*np.pi/180  # deg to rad
    lat2_rad = lat2*np.pi/180
    long2_rad = long2*np.pi/180
    dlat = lat2_rad-lat1_rad
    dlong = long2_rad-long1_rad
    a = (np.sin(dlat/2)**2 +
         np.cos(lat1_rad)*np.cos(lat2_rad)*np.sin(dlong/2)**2)
    c = 2*np.arctan2(np.sqrt(a), np.sqrt(1-a))
    d = R_earth*c  # in km
    return d


def econ_scale_fun(x, C, Q):
    # economy of scale function
    return C * np.power((x/Q), 0.6)


def modified_econ_scale_fun(x, C, Q):
    # modified economy of scale function
    return C + np.power((x/Q), 0.6)


def CapexCurve(xy, fit_type='modified_econ_scale', print_prmtrs=False):
    # returns a capex function. It is calculated with the chosen fit_type
    # using values from the xy dataframe
    x = xy['Capacity_t_y'].values  # t
    y = xy['Capex_USD'].values  # $
    # uses the 4 datapoints in the literature to estimate capex based on capacity
    if fit_type =='linear':  
        prmtrs = np.polyfit(x,y, 1) 
        fun = np.poly1d(prmtrs)
    elif fit_type == '2order':  
        prmtrs = np.polyfit(x, y, 2)
        fun = np.poly1d(prmtrs)
    elif fit_type == 'econ_scale':
        prmtrs, pcov = curve_fit(econ_scale_fun, x, y, bounds=(0, np.inf))
        def fun(x, C=prmtrs[0], Q=prmtrs[1]):
            return C * np.power((x/Q), 0.6)
    elif fit_type == 'modified_econ_scale':
        prmtrs, pcov = curve_fit(modified_econ_scale_fun, x, y, bounds=(0, np.inf))
        def fun(x, C=prmtrs[0], Q=prmtrs[1]):
            return C + np.power((x/Q), 0.6)
    if print_prmtrs:  # return also the prmtrs
        return fun, prmtrs
    else:  # only return the function
        return fun


def FoodWasteDistribution(Spatialdata, feed_names=['SRU', 'DCW', 'BSG'], 
                          Q_min=1.81437, Q_max=40, factor=250,
                          fig_name='FoodDistrMap',
                          subfolder='', transparency=False, ref_40km='circle',
                          plant_LatLonCap=False, size_plant=2400,
                          negative=False,
                          xTicks=[-80, -78, -76, -74, -72], 
                          yTicks=[41, 43, 45],
                          cities = ['White\nPlains',
                                    'Buffalo', 'Syracuse', 'Rochester', 'Albany'], 
                          lat_cities = [41.9, 42.5, 43.4, 42.8, 43.4], 
                          long_cities = [-74.7, -78.7, -75.7, -77.6, -74.3],
                          legend=True):
    """
    Creates a spatial map of NYS food waste producers, differentiated by type
    (feed names) and flowrate.
    the map can include HTC plant locations and capacity if provided

    Parameters
    ----------
    Spatialdata : pandas dataframe data must be given in t/week
        DESCRIPTION.
    feed_names : list of str
        DESCRIPTION.
    Q_min : float, optional
        minimum FW flowrate that is displayed in the map (t/week).
    Q_max : TYPE, optional
        producers with Q>Q_max are represented on the map with a dot of the 
        size of Q_max. (t/week)
    factor : float, optional
        ratio to scale the dot size on the map. The default is 250.
    fig_name : str, optional
        plot file name. The default is 'FoodDistrMap'.
    subfolder : str, optional
        name of subfolder where outputs are saved. The default is ''.
    transparency : bool, optional
        to save png with or without transparency. The default is False.
        ref_40km: bool. If True print a segment that is 40 km long on the map
    plant_locs : list of tuples with latlong data of plants to be plotted
    negative : used to create negative transparent figures
    cities, lat_cities, long_cities: data used to add city names in specific
          locations, 

    """
    
    
    data_Q = [pd.DataFrame()]*len(feed_names)
    fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=0)
    clr = 'k'  
    if negative:  # default text color becomes white
        clr = 'w'
        fig_name = fig_name + '_neg'
        
        transparency = True
    # create a map of NYS
    ax[0] = NY.boundary.plot(color=clr, figsize=(8, 8), aspect='1.4')
    if negative:
        ax[0].set_axis_off()  # do not put x and y axis
    for f, name in enumerate(feed_names):  # filter and add FW producers
        data_Q[f] = Spatialdata[f].loc[Spatialdata[f]['EstWasteTonWeek']
                                              >= Q_min].copy()
        corr_Q = np.where(data_Q[f]['EstWasteTonWeek'] < Q_max,
                          data_Q[f]['EstWasteTonWeek'], Q_max)        
        s_factor = corr_Q/np.max(corr_Q)*factor
        ax[0].scatter(data_Q[f].Longitude, data_Q[f].Latitude,
                      color=clrs[f], s=s_factor, alpha=0.3)
        ax[0].scatter(data_Q[f].Longitude, data_Q[f].Latitude, s=s_factor, 
                      alpha=1, facecolor='none', edgecolor=clrs[f])

    # create entries for the legend    
    SRU = plt.scatter([],[], s=np.max(s_factor)/4*3, color=clrs[0])
    DCW = plt.scatter([],[], s=np.max(s_factor)/4*3, color=clrs[1])
    BSG = plt.scatter([],[], s=np.max(s_factor)/4*3, color=clrs[2])
    mrk0 = plt.scatter([],[], s=np.max(s_factor)/2, color='None')
    mrk1 = plt.scatter([],[], s=np.min(s_factor), color='grey')
    mrk2 = plt.scatter([],[], s=np.max(s_factor)/4, color='grey')
    mrk3 = plt.scatter([],[], s=np.max(s_factor)/2, color='grey')
    mrk4 = plt.scatter([],[], s=np.max(s_factor)/4*3, color='grey')
    mrk5 = plt.scatter([],[], s=np.max(s_factor), color='grey')
    if not plant_LatLonCap:  # only FW map is plotted
        leg_entries = [SRU, DCW, BSG, mrk0, mrk0, mrk1, mrk2, mrk3, mrk4, mrk5]   
        leg_lbls = feed_names + \
            ['', '', str(np.round(Q_min, 1)) + ' t/week', 
              str(int(np.max(corr_Q)/4)), str(int(np.max(corr_Q)/2)),
              str(int(np.max(corr_Q)/4*3)), str(np.round(Q_max, 1)) + '+']
        leg_loc = 'lower left'
        leg_ncol = 2
        xLim=None        
    else:  # best plant locations are provided
        lat_plants = np.asarray([plant[0] for plant in plant_LatLonCap])
        long_plants = np.asarray([plant[1] for plant in plant_LatLonCap])
        cap_plants = np.asarray([plant[2] for plant in plant_LatLonCap], dtype=int)  # t/y
        st_cap_plants = [str(int(c)) for c in cap_plants]  # t/day (Q_tot)
        
        ax[0].scatter(long_plants, lat_plants, s=size_plant, facecolor='none',
                      edgecolor=clrs[3])
        ax[0].plot(long_plants, lat_plants, marker='*', markersize=10,
                   linestyle='None', color=clrs[3])
        # for p in range(len(st_cap_plants)):
        #     ax[0].text(long_plants[p], lat_plants[p] - .1, st_cap_plants[p], 
        #                color=clrs[3],
        #                va='top', ha='center', size=12, weight='bold')
        plnts = plt.scatter([],[], marker='*', s=50, linestyle='None',
                            color=clrs[3])
        leg_entries = [SRU, DCW, BSG, mrk0, mrk1, mrk2, mrk3, mrk4, mrk5,
                       mrk0, plnts]   
        leg_lbls = feed_names + \
            ['', str(np.round(Q_min, 1)) + ' t/week', 
              str(int(np.max(corr_Q)/4)), str(int(np.max(corr_Q)/2)),
              str(int(np.max(corr_Q)/4*3)), str(np.round(Q_max, 1)) + '+',
              '', 'HTC plant\n(capacity \n[t/day])']
        leg_loc = 'upper right'
        leg_ncol = 1
        xLim=[-79.8, -71.45]
    if legend:  # put legend and reference for 40.2 km (either simple or circle) 
        ax[0].legend(handles=leg_entries, labels=leg_lbls, loc=leg_loc,
                     scatterpoints=1, ncol=leg_ncol)
        if ref_40km == 'simple':  
            ax[0].plot([-71.75 - 0.47664, -71.75], [40.5, 40.5], marker='|',
                       linestyle=lnstls[0], color='k')        
            ax[0].text(-71.75 - 0.47664/2, 40.5, '40 km', va='bottom',
                       ha='center', size=10)
        elif ref_40km == 'circle':    
            x, y = -72.3, 41.7
            ax[0].plot([x, x + 0.47664], [y, y], marker='|', color=clr)
            ax[0].text(x + 0.47664/2, y, '40.2 km', va='bottom', ha='center', 
                       color=clr, size=9)
            if plant_LatLonCap:
                size_plant -= 200  # the size of the circle changes based on xLim
            ax[0].scatter(x, y, s=size_plant, alpha=1, facecolor='none',
                          edgecolor=clr)
    if cities:        
        for c, city in enumerate(cities):
            ax[0].text(long_cities[c]+.08, lat_cities[c]-.05, city, va='center',
                        ha='center', size=12, weight='bold', color=clr)
    FigSave(fig_name, out_path, fig, ax, axt, fig_par,
            subfolder=subfolder, xLab='Longitude', yLab='Latitude', xLim=xLim,
            xTicks=xTicks, yTicks=yTicks,
            tight_layout=True, legend=None, transparency=transparency)
    

def SingleLoc(subfolder='SingleLocation',
        # input numerical data, the default values are for the base case
        latlong=(42.82828282828283, -73.77777777777777),  # Buffalo
        Q_min=1.81437, buf_radius=40.2335, T_HTC=220, t_HTC=1, T_t_HTC=None,
        Spatialdata=None, Feeds_props=None, HC_yields_all=None,
        HC_HHVs_all=None, Gas_yields_all=None, Utility_factors_all=None,
        CapexFunc=None, feed_names=['SRU', 'DCW', 'BSG'], 
        # constants or fixed parameters
        transp_price_h=100, transp_price_d=2.50, price_HC_USD_t=16.5,
        price_elect = 0.055,  # [$/kWh]
        labor_rate = 27.5, # $/hr avg for chemical plant operators in NYS
        P_water_treatment = 0.86, # [$/m3] based on syracuse water authority
        tipping_fee=78.1, T_mix=25, prjct_lftm=30, cnstrn_prd=3, 
        dscnt_rt=0.07, price_NatGas=6.73/1000,
        truck_liquid_tons=9000,  # 9000 gal/truck->t/truck (m3=t)
        report=False, HX_plt_T=False, plt_cum_NPV=False,
        cost_emiss_brkdwn=False, xLim_cost_emiss=[None, None],
        xTicks_cost_emiss=[None, None],
        file_name=''):
      
    Q_feeds = np.zeros(len(feed_names))
    d_haul = np.zeros(len(feed_names))
    h_haul = np.zeros(len(feed_names))
    cost_transp = 0
    SpDt = Spatialdata  # rename for conciseness
    SpDt_Q = [pd.DataFrame()]*len(feed_names)
    SpDt_Q_r = [pd.DataFrame()]*len(feed_names)
    lat = latlong[0]
    long = latlong[1]
    if T_t_HTC:  # if HTC conditions are given as a single variable
       T_HTC = T_t_HTC[0]
       t_HTC = T_t_HTC[1]
    for f, name in enumerate(feed_names):
        SpDt_Q[f] = SpDt[f].loc[SpDt[f]['EstWasteTonWeek'] >= Q_min].copy()
        dist = latlongdist(lat, long, SpDt_Q[f]['Latitude'],
                           SpDt_Q[f]['Longitude'])
        dist_filt = dist[dist <= buf_radius]
        SpDt_Q_r[f] = SpDt_Q[f].loc[(dist <= buf_radius)]

        Q_feeds[f] = SpDt_Q_r[f]['EstWasteTonWeek'].sum()/7  # ton/week -> ton/day
        # transportation cost estimation
        if name != 'DCW':
            d_haul[f] = np.sum(SpDt_Q_r[f]['EstWasteTonWeek'] * dist_filt/20*52)   # km/y
            cost_transp += d_haul[f]*transp_price_d  # $/y
        else:
            # trip/week  # gal/truck->t/truck (m3=t)
            trip_week = SpDt_Q_r[f]['EstWasteTonWeek']/truck_liquid_tons*cf_gal_m3
            # [h] truck_vol_gal/pump_gpm/min_hours/load_unload + dist/speed
            hour_trip = (9000/1300/60)*2 + dist_filt/80
            h_haul[f] = np.sum(hour_trip*trip_week)*52  # hour/y
            d_haul[f] = np.sum(trip_week*dist_filt)*52  # km/y
            cost_transp += h_haul[f]*transp_price_h  # $/y
    Q_tot = np.sum(Q_feeds)  # ton/day
    # list of names of variables to store in the output dataframe
    res_labels = \
        ['long', 'lat', 'Q_min', 'buf_radius', 'T_HTC', 't_HTC',  # info
         'Q_tot',  # results before if (includes Q_feeds)
         'Q_HTC_tot', 'HC_gen', 'Q_AP', 'Gas_gen', 'En_gen', 'ER', 'Q_vol', # energy and mass flows
         'V_reactor', 'PlantCapacity', # reacotr sizes 
         'FC', 'C_labor', 'supervision', 'direct_overhead', 'general_overhead',  # fixed costs
         'mainteneance_labor', 'maintenance_material', 'insurance_tax',
         'VC', 'cost_transp', 'cost_HTC_heat', 'cost_HTC_elect',  # variable costs
         'Capex', 'WorkCapital', 'Opex', 'Revenues_TF', 'Revenues_HC', # DCF
         'Norm_NPV', 'NPV', 'LCOE', 'MinSellP_HC', 'MinSellP_HC_eff',
         'Net_CO2_emiss', 'CI_HC_energy', 'CI_HC_mass', 'Carbon_Price'] # emissions
    # feedstock single flows are added in the dataframe
    [res_labels.insert(6+i, 'Q_feeds_' + str(i+1)) for i in range(len(Q_feeds))]
    if Q_tot == 0:  # location without any FW
        res_vals = [np.nan]*len(res_labels)
        res_vals[0:6] = [long, lat, Q_min, buf_radius, T_HTC, t_HTC]
        cost_emiss_brkdwn= False  # no need to evaluate anything without FW
    else:
        Q_perc = Q_feeds/Q_tot  # [-]

        # selection of feed and HC properties for the selected HTC conditions,
        # create vectors with Feed and HC properties for each location
        HC_yields_h = HC_yields_all.loc[HC_yields_all['hours'] == t_HTC]  # [-]
        HC_HHVs_h = HC_HHVs_all.loc[HC_HHVs_all['hours'] == t_HTC]  # [MJ/kg]
        Gas_yields_h = Gas_yields_all.loc[Gas_yields_all['hours'] == t_HTC] # [-]
        Utility_factors_h = Utility_factors_all.loc[Utility_factors_all['hours'] == t_HTC] # [-]

        HC_yields_hT = [HC_yields_h.loc[HC_yields_h['Feed'] == f][T_HTC]
                        for f in feed_names]  # [-]
        HC_HHVs_hT = [HC_HHVs_h.loc[HC_HHVs_h['Feed'] == f][T_HTC]
                      for f in feed_names]  # [MJ/kg]
        Gas_yields_hT = [Gas_yields_h.loc[Gas_yields_h['Feed'] == f][T_HTC]
                      for f in feed_names]  # [-]
        
        feed_densities = [Feeds_props.loc[Feeds_props['Feed'] == f]['density']
                          for f in feed_names]  # [kg/m3]
        feed_mc = [Feeds_props.loc[Feeds_props['Feed'] == f]['moisture']
                   for f in feed_names]  # [-]
        Heat_factor_hT = Utility_factors_h.loc[Utility_factors_h['Utility'] 
                                               == 'Heat'][T_HTC].values[0]  # kWh/kg feed
        Elec_factor_hT = Utility_factors_h.loc[Utility_factors_h['Utility'] 
                                               == 'Electricity'][T_HTC].values[0]  # kWh/kg feed
        # =============================================================================
        # # computation of mixture properties based on feed streams and properties
        # =============================================================================
        mix_densities = np.average(
            feed_densities, weights=Q_perc, axis=0)[0]  # [kg/m3]
        mix_mc = np.average(feed_mc, weights=Q_perc, axis=0)[0]  # [-]
        mix_HC_HHVs = np.average(HC_HHVs_hT, weights=Q_perc, axis=0)[0] # [MJ/kg]
        mix_HC_yields = np.average(HC_yields_hT, weights=Q_perc, axis=0)[0] # [-]
        mix_Gas_yields = np.average(Gas_yields_hT, weights=Q_perc, axis=0)[0]  # [-]
        mix_target_mc = np.where(mix_mc > .85, mix_mc, .85)  # [-]
        mix_water_add = (Q_tot*(mix_target_mc - mix_mc) / 
                         (1 - mix_target_mc)) # [t/day]
        # mix_water_tot = mix_water_add + Q_tot*mix_mc
        # BWratio = mix_dry_matter/mix_water_tot  # dry biomass/water ratio
        # total flow entering HTC reactor [t/day]
        # Q_tot_dry = Q_tot*(1 - mix_mc)  # [t/d] dry matter in Q_tot
        Q_HTC_tot = Q_tot + mix_water_add  # [t/d]
        HC_gen = mix_HC_yields*Q_HTC_tot*(1-mix_target_mc)  # HC generation [t/day]
        En_gen = HC_gen*mix_HC_HHVs*1e3  # total energy generation [MJ/d]
        Gas_gen = mix_Gas_yields*Q_HTC_tot * (1-mix_target_mc) # gas generation [t/day]
        # energy recovery HC/feed [-]
        ER = En_gen/(Q_tot*(1-mix_mc)*mix_HC_HHVs*1000)
        Q_vol = Q_HTC_tot*(1/mix_densities)*1000  # volumetric flow [m3/day]
        V_reactor = Q_vol/24*t_HTC  # continuous reactor, filling = 1 [m3]
        Q_AP = Q_HTC_tot - HC_gen - Gas_gen # aquous phase genration (t/day)
        # =====================================================================
        # cost estimation discounted cash flow calculations
        # =====================================================================
        # one time costs
        
        PlantCapacity = Q_HTC_tot*365  # wet t/y used for plant sizing
        Capex = CapexFunc(PlantCapacity)  # [$] uses total dry basis
        WorkCapital = 0.05*Capex  # working capital  # [$]  
        # annual cost and labor computanion
        HTC_production = (HC_gen*1000/24)  # (kg/h)        
        # annual HTL labor hours. 1 unit process (info in SI)
        HTC_labor_hours = 2.13*(HTC_production**0.242)*1*(8760/24) # [h/y]
        C_labor = HTC_labor_hours*labor_rate  # [$/y]
        supervision = 0.25*C_labor
        direct_overhead = 0.5*(C_labor+supervision)
        general_overhead = 0.5*(C_labor+supervision+direct_overhead)
        mainteneance_labor = 0.015*Capex
        maintenance_material = 0.015*Capex
        insurance_tax = 0.01*Capex
        FC = (C_labor + supervision + direct_overhead + general_overhead +
              mainteneance_labor + maintenance_material + insurance_tax)
        # HTC_heat [t/d]*[kg/t]*[d/y]*[kWh/kg]*[MJ/kWh]=[MJ/y] (info SI)
        HTC_heat = (Q_HTC_tot*1000)*365*Heat_factor_hT*3.6  
        # price_NG_MJ = price per ft3*conv_ft_m3/density_CH4/HHV_CH4
        price_NG_MJ = price_NatGas*35.3147/0.668/55  # [$/MJ]
        cost_HTC_heat = HTC_heat*price_NG_MJ
        # from fiori paper
        
        # HTC_elect_consumpt [t/d]*[kg/t]*[d/y]*[kWh/kg]=[kWh/y] (info SI)
        HTC_elect_consumpt = (Q_HTC_tot*1000)*365*Elec_factor_hT  # [kWh/y]
        cost_HTC_elect = HTC_elect_consumpt * price_elect  # [$/y]
        
        # Aqueous phase treatment costs
        # cost_AP_treatment [t/day]*[kg/t]*[m3/kg]*[$/m3]*[day/y] = [$/y]
        cost_AP_treatment = (Q_AP*1000/1000)*P_water_treatment*365 
        VC = (cost_transp + cost_HTC_heat + cost_HTC_elect
              + cost_AP_treatment)  # [$/y]
        Opex = VC + FC  # [$/y]
        # =====================================================================
        # DCF computation  
        # =====================================================================
        Q_tot_y = Q_tot*365  # t/day --> t/y
        HC_gen_y = HC_gen*365  # t/day --> t/y
        En_gen_y = En_gen*365  # [MJ/day]*day/y --> [MJ/y]
        Revenues_TF = Q_tot_y*tipping_fee  # t/y * $/t -> $/y
        Revenues_HC = HC_gen_y*price_HC_USD_t  # -> $/y
        
        prjct_lftm += 1 # project lifetime - add 1 due to python 0 to x notation
        # actualization factor
        actlztn_fctr = (1 + dscnt_rt)**np.arange(prjct_lftm)  # [-]

        CAPEX = np.zeros(prjct_lftm)  #  [$]
        OPEX = np.zeros(prjct_lftm)  # [$/y]
        REVENUES_TF = np.zeros(prjct_lftm)  # [$/y]
        REVENUES_HC = np.zeros(prjct_lftm)
        WC = np.zeros(prjct_lftm)  # [$]
        
        FIXED_COSTS = np.zeros(prjct_lftm)
        COST_TRANS = np.zeros(prjct_lftm)
        COST_HTC_HEAT= np.zeros(prjct_lftm)
        COST_HTC_ELECT = np.zeros(prjct_lftm)
        COST_AP_TREATMENT = np.zeros(prjct_lftm)
        FCF = np.zeros(prjct_lftm)  # [$/y]
        PV = np.zeros(prjct_lftm)  # [$]
        Cost_discounted = np.zeros(prjct_lftm)   # [$/y]
        # Q_HTC_discounted = np.zeros(prjct_lftm)  # t/y
        En_gen_y_discounted = np.zeros(prjct_lftm)  # [MJ/y]

        
        CAPEX[0] = Capex  # [$] 
        WC[cnstrn_prd] = WorkCapital  # [$]
        OPEX[cnstrn_prd+1:] = Opex  # [$/y]        
        REVENUES_TF[cnstrn_prd+1:] = Revenues_TF # [$/y]
        REVENUES_HC[cnstrn_prd+1:] = Revenues_HC # [$/y]
        FIXED_COSTS[cnstrn_prd+1:] = FC # [$/y]
        COST_TRANS[cnstrn_prd+1:] = cost_transp # [$/y]
        COST_HTC_HEAT[cnstrn_prd+1:] = cost_HTC_heat # [$/y]
        COST_HTC_ELECT[cnstrn_prd+1:] = cost_HTC_elect # [$/y]
        COST_AP_TREATMENT[cnstrn_prd+1:] = cost_AP_treatment # [$/y]
        # revenues total
        FCF = (REVENUES_TF + REVENUES_HC) - (CAPEX + WC + OPEX)
        PV = FCF/actlztn_fctr
        # cum_NPV = np.cumsum(PV)
        NPV = np.sum(PV)

        Q_HTC_discounted = Q_tot_y/actlztn_fctr[cnstrn_prd+1:]  # [t/y]
        HC_gen_discounted = HC_gen_y/actlztn_fctr[cnstrn_prd+1:]  # [t/y]
        En_gen_y_discounted = En_gen_y/actlztn_fctr[cnstrn_prd+1:]  # [MJ/y]
        Revenues_TF_discounted = REVENUES_TF/actlztn_fctr  # [$/y]
        
        Norm_NPV = NPV/np.sum(Q_HTC_discounted)  # $/t
        Cost_discounted = (CAPEX + WC + OPEX)/actlztn_fctr
        LCOE = (np.sum(Cost_discounted)/
                np.sum(En_gen_y_discounted))  # $/MJ
        MinSellP_HC = (np.sum(Cost_discounted)/
                       np.sum(HC_gen_discounted))  # $/t
        
        Norm_Revenues_TF = np.sum(Revenues_TF_discounted)/np.sum(HC_gen_discounted)
        MinSellP_HC_eff = MinSellP_HC - Norm_Revenues_TF  # $/t
                
        # =====================================================================
        # GHG emissions
        # =====================================================================
        # Avoided landfill emissions
        EF_landfills = 493/cf_st_t # kg CO2eq/ton -> kg CO2eq/t
        EF_coalmining = 74 # kg CO2eq/t
        CH4_EF_transportation = 0.0051/1000/cf_mile_km # g/mile -> kg/km
        N2O_EF_transportation = 0.0048/1000/cf_mile_km # g/mile -> kg/km
        CO2_EF_transportation = 10.21/cf_gal_m3 # kg/gal diesel -> kg/m3
        fuel_economy = 6.6*cf_mile_km/cf_gal_m3 # mile/gal (heavy duty truck) -> km/m3
        GWP_CH4 = 28.0
        GWP_N20 = 265.0
        EF_NG = 51.66/1000 # kgCO2/MJ
        # (from eGRID database, info in SI)
        EF_electric_NYgrid = 416.69*cf_pound_kg # lb CO2eq/MWh  -> kgCO2eq/MWh
        
        CO2_landfills = Q_tot*EF_landfills*365 # kg CO2eq/y
        # Avoided coal mining emissions (coal displacement with hydrochar)
        
        CO2_coalmining = HC_gen*EF_coalmining*365 # kg CO2eq/y
        # Transportation emissions
        
        CO2_transportation = ((CH4_EF_transportation*GWP_CH4+
                              N2O_EF_transportation*GWP_N20+
                              CO2_EF_transportation/fuel_economy)
                              *np.sum(d_haul)) # kg CO2eq/yr
        # heat and electric import emissions
        
        CO2_heat = EF_NG*HTC_heat # kg CO2/y
        CO2_electric = EF_electric_NYgrid*HTC_elect_consumpt/1000 # kg CO2eq/yr
        CO2_HTC = Gas_gen*365*1000  # t/day -> t/y -> kg/y
        # net GHG emissions and carbon intensity
        Net_CO2_emiss = (CO2_heat+CO2_electric+CO2_transportation+CO2_HTC
                         -CO2_coalmining-CO2_landfills)/1000 # t CO2eq/yr
        CI_HC_energy = (Net_CO2_emiss*1e6)/En_gen_y # g CO2eq/MJ
        CI_HC_mass = (Net_CO2_emiss)/HC_gen_y # t CO2eq/t
        
        # Competitive carbon price 
        if CI_HC_mass < 0:
            carbon_price = MinSellP_HC_eff/abs(CI_HC_mass) # [$/t CO2eq]
        else:  # if CI > 0 there is no carbon saving
            carbon_price = 0

        res_vals = [long, lat, Q_min, buf_radius, T_HTC, t_HTC,  # info
                      Q_tot, # results before if (includes Q_feeds)
                      Q_HTC_tot, HC_gen, Q_AP, Gas_gen, En_gen, ER, Q_vol, # energy and mass flows
                      V_reactor, PlantCapacity, # reactor sizes
                      FC, C_labor, supervision, direct_overhead, general_overhead,  # fixed costs
                      mainteneance_labor, maintenance_material, insurance_tax,
                      VC, cost_transp, cost_HTC_heat, cost_HTC_elect,  # variable costs
                      Capex, WorkCapital, Opex, Revenues_TF, Revenues_HC, # DCF
                      Norm_NPV, NPV, LCOE, MinSellP_HC, MinSellP_HC_eff,
                      Net_CO2_emiss, CI_HC_energy, CI_HC_mass, carbon_price] # emissions
        [res_vals.insert(6+i, Q_feeds[i]) for i in range(len(Q_feeds))]
        res_vals = [float(d) for d in res_vals]

    df_vrbls = pd.DataFrame(res_vals).T

    df_vrbls.columns = res_labels
    if cost_emiss_brkdwn:  # cost and emissions brkdwn
        # checks and creates the folder for outputs if missing
        plib.Path(plib.Path(out_path, subfolder)).mkdir(parents=True,
                                                        exist_ok=True)
        # cost breakdown
        cost_names = ['Tipping Fees (+)', 'Selling HC (+)', 'CAPEX',  
                      'Fixed Costs (-)',  'HTC heat', 'HTC electricity',
                      'Aq. phase treatment', 'Transportation', 
                      'Working Capital']
        costs_not_act = np.asarray([REVENUES_TF, REVENUES_HC, -CAPEX, 
                                    -FIXED_COSTS, -COST_HTC_HEAT,
                                    -COST_HTC_ELECT, -COST_AP_TREATMENT, 
                                    -COST_TRANS, -WC])  # [$]
        cost_acts = np.zeros(len(costs_not_act))
        for v, vec in enumerate(costs_not_act):
            cost_acts[v] = np.sum(vec/actlztn_fctr)
        df_cost_brkdwn = pd.DataFrame(data=dict(zip(cost_names, 
                                                    cost_acts)), index=[0])
          
        # emission breakdow
        v_name_emiss = ['Landfills', 'Coal Mining', 'Heat', 'Electricity', 
                        'HTC gas', 'Transportation' ]
        v_emiss = np.asarray([-CO2_landfills, -CO2_coalmining,
                              CO2_heat, CO2_electric,
                              CO2_HTC, CO2_transportation]) # kg/y
        
        df_emiss_brkdwn = pd.DataFrame(data=dict(zip(v_name_emiss, 
                                                     v_emiss)), index=[0])

        
        fig_name = 'CostEmiss_' + file_name          
        y_costs= np.linspace(len(cost_acts)-1, 0, len(cost_acts))
        y_emiss= np.linspace(len(v_emiss)-1, 0, len(v_emiss))
        fig, ax, axt, fig_par = FigCreate(rows=2, cols=1, plot_type=0,
                                          paper_col=.78, hgt_mltp=1.25)
        ax[0].barh(y_emiss, v_emiss/1e6, color=clrs, edgecolor='k') 
        for e in range(len(v_emiss)):
            if v_emiss[e] > 0:
                ax[0].text(-1.2, y_emiss[e], v_name_emiss[e], va='center',
                            ha='right')
            else:
                ax[0].text(1.2, y_emiss[e], v_name_emiss[e], va='center',
                            ha='left')
                
        ax[1].barh(y_costs, cost_acts/1e6, color=clrs, edgecolor='k') 
        for c in range(len(cost_acts)):
            if cost_acts[c] > 0:
                ax[1].text(-2, y_costs[c], cost_names[c], va='center',
                           ha='right')
            else:
                ax[1].text(2, y_costs[c], cost_names[c], va='center',
                           ha='left')

        FigSave(fig_name, out_path, fig, ax, axt, fig_par,
                subfolder=subfolder,
                xLim=xLim_cost_emiss, xTicks=xTicks_cost_emiss,
                yTicks=[[], []], yTickLabels=[[], []],
                xLab=['Emissions [kt CO$_{2, eq.}$/y]', 'Revenues/Costs [M$]'],
                legend=None, annotate_lttrs=True, tight_layout=True)
    if report:
        df_cost_brkdwn.to_csv(plib.Path(out_path, subfolder,
                                        'CostsUSD' + file_name + '.csv')) 
        df_emiss_brkdwn.to_csv(plib.Path(out_path, subfolder, 
                                         'EmisskgCO2_y' + file_name + '.csv'))
    return df_vrbls


def MultiLoc(latslongs, prmtrs, vrbls, name='name',
             feed_names=['SRU', 'DCW', 'BSG'],
             plt_cum_NPV=False, print_report=True, df_to_csv=True,
             plt_grid=True, plt_grid_ax_latlong=True, print_prmtrs=True,
             subfolder='HeatMap', best_plt=False):
    combs = list(itr.product(latslongs, prmtrs[0].value, prmtrs[1].value,
                             prmtrs[2].value, prmtrs[3].value))
    lst_dfs = []
    for comb in combs:
        lst_dfs.append(SingleLoc(latlong=comb[0], Q_min=comb[1],
                                 buf_radius=comb[2], T_HTC=comb[3],
                                 t_HTC=comb[4],
                                 feed_names=feed_names,
                                 Spatialdata=Spatialdata,
                                 Feeds_props=Feeds_props,
                                 HC_yields_all=HC_yields_all,
                                 HC_HHVs_all=HC_HHVs_all,
                                 Gas_yields_all=Gas_yields_all,
                                 Utility_factors_all=Utility_factors_all,
                                 CapexFunc=CapexFunc,
                                 report=False, plt_cum_NPV=plt_cum_NPV))
    df = pd.concat(lst_dfs, ignore_index=True)
    if not best_plt:
        if prmtrs and all([pr.type == 'number' for pr in prmtrs]):
            name = 'HeatMap_' + '_'.join([str(prmtr.value[0])
                                          for prmtr in prmtrs])
        if df_to_csv:  # save files to csv
            plib.Path(plib.Path(out_path, subfolder)).mkdir(parents=True,
                                                            exist_ok=True)
            csv_path = plib.Path(out_path, subfolder, name + '.csv')
            df.to_csv(csv_path)
        if plt_grid:
            GridPlot(df, vrbls=vrbls, prmtrs=prmtrs, subfolder=subfolder,
                     name=name, best_plt=best_plt, 
                     ax_latlong=plt_grid_ax_latlong, print_prmtrs=print_prmtrs)
    else:
        if df_to_csv:  # save files to csv
            plib.Path(plib.Path(out_path, subfolder, name)).mkdir(parents=True, 
                                                                  exist_ok=True)
            csv_path = plib.Path(out_path, subfolder, name, name + '.csv')
            df.to_csv(csv_path)
        if plt_grid:
            GridPlot(df, vrbls=vrbls, prmtrs=prmtrs, subfolder=subfolder,
                     name=name, best_plt=best_plt,
                     ax_latlong=plt_grid_ax_latlong, print_prmtrs=print_prmtrs)
    return df


def BestLocSelect(grid_df, prmtrs, vrbls, name='name', choice_var='NPV',
                  choice_maximize=True,
                  n_best=5, skip_n_loc=1, min_dist_best=40.2335*2, df_best_to_csv=True,
                  plot=True, plt_grid_ax_latlong=True, subfolder='BL',
                  transparency=False, negative=False,
                  cities = ['NYC', 'Buffalo', 'Syracuse',
                            'Rochester', 'Albany'],
                  lat_cities = [40.9, 42.9, 43.4, 43.5, 43.35],
                  long_cities = [-74.6, -79.5, -75.7, -77.4, -74.25],
                  legend=True):
    max_computation = np.max(grid_df.index)
    if choice_maximize:  # order from max var to min var
        sort_grid_df = grid_df.sort_values(by=[choice_var], ascending=False,
                                           ignore_index=True)
    else:  # order from min var to max var
        sort_grid_df = grid_df.sort_values(by=[choice_var], ignore_index=True)
    # becomes the list of indexes of best locs
    df_best = sort_grid_df.iloc[[0]]
    c = 1
    while (len(df_best) < n_best + skip_n_loc) and (c < max_computation):
        # go thru the best locations until the next non-close one is found
        condition = True
        while condition:
            dist = latlongdist(sort_grid_df.lat.iloc[c],
                               sort_grid_df.long.iloc[c],
                               df_best.lat, df_best.long)
            
            if all(dist >= min_dist_best):
                condition = False
                df_best = pd.concat([df_best, sort_grid_df.iloc[[c]]],
                                    ignore_index=True)
            c += 1
    df_best = df_best.iloc[skip_n_loc:]
    df_best = df_best.reset_index(drop=True)
    if prmtrs and all([prmtr.type == 'number' for prmtr in prmtrs]):
        name = 'BestLoc_' + '_'.join([str(prmtr.value[0])
                                      for prmtr in prmtrs]) + '_' + choice_var
    if df_best_to_csv:  # save files to csv
        # plib.Path(plib.Path(out_path, subfolder, name)).mkdir(parents=True, 
                                                              # exist_ok=True)
        csv_best_path = plib.Path(out_path, subfolder, name + '.csv')
        df_best.to_csv(csv_best_path)
    if plot:
        plant_LatLonCap = []
        for p in range(df_best.index[-1] + 1):
            plant_LatLonCap.append((df_best.lat[p], 
                                    df_best.long[p],
                                    df_best.Q_tot[p]))
        FoodWasteDistribution(Spatialdata, feed_names,
                              plant_LatLonCap=plant_LatLonCap,
                              subfolder=subfolder,transparency=transparency,
                              negative=negative,
                              fig_name='FoodDitrMap_' + name,
                              cities=cities, lat_cities=lat_cities,
                              long_cities=long_cities, legend=legend)
    return df_best


def GridPlot(res, vrbls, prmtrs=None, subfolder='',
              name='name', choice_var=None, best_plt=False,
              ax_latlong=True, cmap='coolwarm', print_prmtrs=True,
              annotate_lttr=False):    
    if choice_var or best_plt:  # means only best locations are plotted
        s_size = 50  # max(res.index()    
        subfolder = subfolder + '\\' + name        
    else:  # heatmap
        s_size = 10
    for vrbl in vrbls:
        print('Plotting', vrbl.name)        
        data = res[vrbl.name].values*vrbl.conv_factor
        if vrbl.limit:
            data = np.where(data < vrbl.limit[0], vrbl.limit[0], data)
            data = np.where(data > vrbl.limit[1], vrbl.limit[1], data)
        fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=0)
        ax[0] = NY.boundary.plot(color='k', figsize=(8, 8), aspect='1.4')
        cb = ax[0].scatter(res.long, res.lat, c=data, s=s_size, cmap=cmap)        
        if print_prmtrs and all([pr.type == 'number' for pr in prmtrs]):
            text = ''.join([prmtr.label + ' = ' + 
                            str(np.round(prmtr.value[0], 1)) +
                            ' ' + prmtr.unit + '\n ' for prmtr in prmtrs])
            ax[0].text(-79.6, 41, text, va='center', ha='left', size=13)
        if ax_latlong:
            xLab='Longitude'
            yLab='Latitude'
            ax[0].set_xticks([-80, -78, -76, -74, -72])
            ax[0].set_yticks([41, 43, 45])
            cbar = plt.colorbar(cb, shrink=.6, anchor=(0, .7), pad=-.17)
            cbar.set_label(label=vrbl.label)

        else:
            xLab=None
            yLab=None
            ax[0].spines["top"].set_visible(False)
            ax[0].spines["right"].set_visible(False)
            ax[0].spines["bottom"].set_visible(False)
            ax[0].spines["left"].set_visible(False)
            ax[0].set_yticks([])
            ax[0].set_xticks([])
            cbar = plt.colorbar(cb, shrink=.6, anchor=(0, .7), pad=-.17)
            cbar.set_label(label=vrbl.label)
        if vrbl.limit:  # put < or > at first and last ticklabel
            cbar.ax.set_ylim([vrbl.limit[0], vrbl.limit[1]])
            cbar_ticks = cbar.get_ticks()
            cbar.set_ticks(cbar_ticks)
            nice_ticks = [str(int(c)) for c in cbar_ticks]
            if vrbl.pm_sign[0]:
                nice_ticks[-1] = '> ' + nice_ticks[-1]
            if vrbl.pm_sign[1]: 
                nice_ticks[0] = '< ' + nice_ticks[0]
            cbar.set_ticklabels(nice_ticks)
        if vrbl.annotate_lttr:
            ax[0].annotate('(' + vrbl.annotate_lttr + ')',
                          xycoords='axes fraction',
                          xy=(0, 0), rotation=0, size='x-large',
                          xytext=(0, 0.02), weight='bold')
        FigSave(name + '_' + vrbl.name, out_path, fig, ax, axt, fig_par,
                subfolder=subfolder, 
                xLab=xLab, yLab=yLab,
                legend=None, tight_layout=True) 

#%%
     
def CumulSystemPerf(grid_df, prmtrs, vrbls, name='name', choice_var='NPV',
                    choice_maximize=True, n_best=13, min_dist_best=40.2335*2,
                    skip_n_loc=1,
                    df_best_to_csv=True, plot=True, 
                    subfolder='CumulSystemPerf', file_name='CSP',
                    yLim=None, ytLim=None,
                    xTicks=None, yTicks=None,
                    ytTicks=None, 
                    legend='lower center',
                    clr_axis_labels=False, axis_labels=True, paper_col=.8):
    fig_name = [file_name] 
    for vrb in vrbls:
        fig_name.append(vrb.name)
        
    max_computation = np.max(grid_df.index)
    if choice_maximize:  # order from max var to min var
        sort_grid_df = grid_df.sort_values(by=[choice_var], ascending=False,
                                           ignore_index=True)
    else:  # order from min var to max var
        sort_grid_df = grid_df.sort_values(by=[choice_var], ignore_index=True)
    # becomes the list of indexes of best locs
    df_best = sort_grid_df.iloc[[0]]
    c = 1
    while (len(df_best) <= n_best + skip_n_loc) and (c < max_computation):
        # go thru the best locations until the next non-close one is found
        condition = True
        while condition:
            dist = latlongdist(sort_grid_df.lat.iloc[c],
                               sort_grid_df.long.iloc[c],
                               df_best.lat, df_best.long)
            if all(dist >= min_dist_best):
                condition = False
                df_best = pd.concat([df_best, sort_grid_df.iloc[[c]]],
                                    ignore_index=True)
            c += 1
    df_best = df_best.iloc[skip_n_loc:]
    df_trends = pd.DataFrame(columns=list(grid_df))
    n_vec = np.arange(1, n_best + 1)
    for n in n_vec:
        df_trends.loc[str(n)] = df_best.head(n).sum(numeric_only=True, axis=0)
    fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=1,
                                      paper_col=paper_col, hgt_mltp=1.25)  
    ax[0].plot(n_vec, df_trends[vrbls[0].name]*vrbls[0].conv_factor,
               color=clrs[0], marker=mrkrs[0], linestyle='none',
               label=vrbls[0].label_no_unit)   
    axt[0].plot(n_vec, df_trends[vrbls[1].name]*vrbls[1].conv_factor, 
                color=clrs[1], marker=mrkrs[1], linestyle='none', 
                label=vrbls[1].label_no_unit)  
    if axis_labels:
        xLab='# of HTC plants in NY' 
        yLab=vrbls[0].label
        ytLab=vrbls[1].label
        if clr_axis_labels:
            ax[0].yaxis.label.set_color(clrs[0])
            axt[0].yaxis.label.set_color(clrs[1])
    else:
        xLab, yLab, ytLab = None, None, None
    
    FigSave('_'.join(fig_name), out_path, fig, ax, axt, fig_par,
            subfolder=subfolder, 
            xLim=[1, n_best], yLim=yLim, ytLim=ytLim,
            xTicks=xTicks, yTicks=yTicks, ytTicks=ytTicks,
            xLab=xLab, yLab=yLab, ytLab=ytLab,
            legend=legend, tight_layout=True) 

    return df_trends, df_best

       
def SensAnalysis(prmtrs, vrbls, subfolder='SensitivityAnalysis',
                 xLab='Parameters levels [-]', paper_col=1,
                 file_name='',
                 type_analysis='multi_param', 
                 plt_delta=True, xLim=None, plt_prmt=None,
                 yLim=None, ytLim=None, yTicks=None, ytTicks=None,
                 legend='best', ncol_leg=2):
    # variable (vrbl) are the final value we are interested in
    # parameter (pr) are the ones changed to assess the effect on the vars
    # value (vl)
    # check folder existence and possibly create it if missing
    plib.Path(plib.Path(out_path, subfolder)).mkdir(parents=True, exist_ok=True)
    dfs, dfs_dlt, yLabs, varsLeg,= [], [], [], [] # for iteration over more target variables
    for sens_vrbl in vrbls:  # computation is repeated for each target var
        
        prs_vary = []
        for pr in prmtrs:
            if not pr.do_not_vary:
                prs_vary.append(pr)
        pr_vary_names = [pr.name for pr in prs_vary]
        pr_names = [pr.name for pr in prmtrs]
        pr_labels = [pr.label for pr in prmtrs]
        pr_base_vls = [pr.base_value for pr in prmtrs]
        dc_pr_base = dict(zip(pr_names, pr_base_vls))
        dc_labels = dict(zip(pr_names, pr_labels)) 
        if type_analysis=='multi_param':
            n_steps = prmtrs[0].n_steps  # it is taken from the 1st param
            step_list = np.arange(-(n_steps - 1)/2, (n_steps - 1)/2 + 1,
                                  dtype=int)
            df_pr_steps = pd.DataFrame(index=step_list, columns=pr_names)
            for c, col in enumerate(list(df_pr_steps)):
                df_pr_steps[col] = prmtrs[c].steps
        elif type_analysis=='single_param':
            n_steps = prs_vary[0].n_steps
            step_list = prs_vary[0].steps       
        df_res = pd.DataFrame(index=step_list, columns=pr_vary_names)
        base_case_res = False # to enter the if clause down
        for p, pr in enumerate(prs_vary):  # only calculate pr that vary
            dc_pr_vls = dc_pr_base.copy()  # but use the values of all prs
            for s in range(n_steps):
                dc_pr_vls[pr.name] = pr.steps[s]  # change the pr value
                            
                SL = SingleLoc(\
                       feed_names=feed_names,
                       Spatialdata=Spatialdata,
                       Feeds_props=Feeds_props,
                       HC_yields_all=HC_yields_all,
                       HC_HHVs_all=HC_HHVs_all,
                       Gas_yields_all=Gas_yields_all,
                       Utility_factors_all=Utility_factors_all,
                       CapexFunc=CapexFunc,
                       latlong=dc_pr_vls['latlong'], 
                       T_t_HTC=dc_pr_vls['T_t_HTC'],
                       Q_min=dc_pr_vls['Q_min'],
                       buf_radius=dc_pr_vls['buf_radius'],
                       transp_price_h=dc_pr_vls['transp_price_h'],
                       transp_price_d=dc_pr_vls['transp_price_d'],
                       price_HC_USD_t=dc_pr_vls['price_HC_USD_t'],
                       price_elect=dc_pr_vls['price_elect'],
                       price_NatGas=dc_pr_vls['price_NatGas'],
                       labor_rate=dc_pr_vls['labor_rate'],
                       P_water_treatment=dc_pr_vls['P_water_treatment'],
                       tipping_fee=dc_pr_vls['tipping_fee'],
                       dscnt_rt=dc_pr_vls['dscnt_rt'],
                       report=False, plt_cum_NPV=False,
                       cost_emiss_brkdwn=False)
                df_res.iloc[s][pr.name] = (SL[sens_vrbl.name][0]
                                           *sens_vrbl.conv_factor)
                # base case resutls are stored to compute variations
                if not base_case_res and s==int((n_steps - 1)/2):  
                    
                    base_case_res = (SL[sens_vrbl.name][0]
                                      *sens_vrbl.conv_factor)  
        
        df_var = df_res - base_case_res  # unsorted
        if sens_vrbl.to_maximize:  
            df_var = df_var.T.sort_values(df_var.index[-1], 
                                          ascending=False).T # sorted by impact
        else:
            df_var = df_var.T.sort_values(df_var.index[-1], 
                                             ascending=True).T # sorted by impact 
        df_var.fillna(0)
        dfs.append(df_res)
        dfs_dlt.append(df_var)
        yLabs.append(sens_vrbl.label)
        varsLeg.append(sens_vrbl.label_no_unit)
    
    if type_analysis=='multi_param':
        csv_sens_path = plib.Path(out_path, subfolder, file_name + '.csv')
        df_pr_steps.to_csv(csv_sens_path)
        fig, ax, axt, fig_par = FigCreate(rows=len(dfs_dlt), cols=1, plot_type=0,
                                          paper_col=.95, hgt_mltp=1.5)  
        for d, df in enumerate(dfs_dlt):
            ax[d].hlines(y=0, xmin=-n_steps, xmax=n_steps, color='grey',
                         linestyle=lnstls[1])  # print zero line
            c = 0  # to start line properties always from zero
            for col in list(df):
                if any(df[col] != 0):  # skip prs that have no impact 
                    ax[d].plot(df.index.values, df[col], color=clrs[c], 
                               marker=mrkrs[c], linestyle=lnstls[c],
                               label=dc_labels[col])
                    c += 1  # go to next properties (color, linestyle, mrkr)
        step_list = np.asarray(step_list, dtype=int).tolist()  # avoids problems on MAC OS# 
        FigSave(file_name, out_path, fig, ax, axt, fig_par,
                xLim=[step_list[0], step_list[-1]], xTicks=step_list,
                yLim=yLim, yTicks=yTicks, subfolder=subfolder, xLab=xLab,
                yLab=yLabs, tight_layout=True, annotate_lttrs=True,
                ncol_leg=ncol_leg, legend=legend)
    elif type_analysis=='single_param':  # max of 2 variables
        fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=1,
                                          paper_col=.78, hgt_mltp=1.5)
        col = pr_vary_names[0]      
        ax[0].plot(prs_vary[0].steps, dfs[0][col], color=clrs[0], 
                   linestyle=lnstls[0], label=varsLeg[0])
        axt[0].plot(prs_vary[0].steps, dfs[1][col], color=clrs[1], 
                   linestyle=lnstls[1], label=varsLeg[1])            
        FigSave(file_name, out_path, fig, ax, axt, fig_par,
                xLim=xLim, xLab=dc_labels[col],
                yLim=yLim, ytLim=ytLim, yTicks=yTicks, 
                ytTicks=yTicks, subfolder=subfolder,
                yLab=yLabs[0], ytLab=yLabs[1], 
                tight_layout=True,
                ncol_leg=ncol_leg, legend=legend)    
    return dfs, dfs_dlt 


class prmtr:  # used to store the properties of parameters 
    def __init__(self, name=None, value=None, label=None, unit=None):
        self.name = name
        if value:
            if type(value) == int or type(value) == float:
                self.value = [value]
                self.type = 'number'
            else:
                self.value = value
                self.type = 'list'
        self.label = label
        self.unit = unit


class vrbl:  # used to store the properties of variables (target)
    def __init__(self, name=None, label=None, limit=None, pm_sign=[True, True],
                 label_no_unit=None, conv_factor=1, to_maximize=True,
                 annotate_lttr=False):
        self.name = name
        self.label = label
        self.label_no_unit = label_no_unit
        self.limit = limit
        self.pm_sign = pm_sign
        self.conv_factor = conv_factor
        self.to_maximize = to_maximize
        self.annotate_lttr = annotate_lttr
        
class sens_prmt:  # prmtrs varied in the SensAnalysis
    def __init__(self, name=None, label=None, unit=None, base_value=None, 
                 steps=None, label_no_unit=None,
                 perc_variations=[-.8, -.5, 0, .5, .8],
                 positive=True, do_not_vary=False, n_steps=None):
        self.name = name
        if type(steps) == np.ndarray:
            steps = steps.tolist()
        if not do_not_vary:
            if not steps:  # if no steps are provided
                n_steps = len(perc_variations)
                if not positive:
                    perc_variations = [-p for p in perc_variations]
                steps = base_value*(1 + np.asarray(perc_variations))
            else:
                n_steps = len(steps)
                base_value = steps[int((n_steps-1)/2)]
        self.steps = steps
        self.base_value = base_value     
        self.n_steps = n_steps     
        self.label = label
        self.unit = unit
        self.positive = positive
        self.do_not_vary = do_not_vary
        self.label_no_unit = label_no_unit

# =============================================================================
# # end of function definition
# =============================================================================
# general prmtr definition for plotting purposes
clrs = sns.color_palette('deep', 13)  # list with colors
sns.set("paper", font_scale=1.25)
sns.set_style("ticks")  # sns.set_style("whitegrid")     
# list with linestyles for plotting
lnstls = [(0, ()),  # solid
          (0, (1, 1)),  # 'densely dotted'
          (0, (5, 1)),  # 'densely dashed'
          (0, (3, 1, 1, 1)),  # 'densely dashdotted'
          (0, (3, 1, 1, 1, 1, 1)),  # 'densely dashdotdotted'
          (0, (5, 5)),  # 'dashed'
          (0, (3, 5, 1, 5)),  # 'dashdotted'
          (0, (1, 5)),  # dotted
          (0, (3, 5, 1, 5, 1, 5)),  # 'dashdotdotted'
          (0, (1, 10)),  # 'loosely dotted'
          (0, (5, 10)),  # 'loosely dashed'
          (0, (3, 10, 1, 10)),  # 'loosely dashdotted'
          (0, (3, 10, 1, 10, 1, 10))]  # 'loosely dashdotdotted'
# list with letters for plotting
lttrs = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm']
# list with markers for plotting
mrkrs = ["o", "v", "X", "s", "p", "^", "P", "<", ">", "*", "d", "1", "2", "3"]
in_path, out_path = PathsCreate(__file__)

# conversion factor (all multiply the left unit to obtain the right one)
cf_st_t = 0.907185  # conversion factor from short ton (US) to metric ton (t)
cf_mile_km = 1.60934  # conversion factorfrom mile to km
cf_gal_m3 = 3.785/1000  # conversion factor from gal to m3
cf_pound_kg = 1/2.2046  # conversion factor from pounds to kg 
# start of computations

feed_names = ['SRU', 'DCW', 'BSG']
# lat/long [deg], EstFoodWasteTonWeek [short ton/week]->[t/week]
Spatialdata = [pd.DataFrame()]*len(feed_names)
for f, name in enumerate(feed_names):
    Spatialdata[f] = pd.read_excel(plib.Path(in_path, 'Spatialdata_HTC.xlsx'),
                                   sheet_name=name, engine='openpyxl')
    Spatialdata[f].EstWasteTonWeek *= cf_st_t  
# Feed_props: HHV [MJ/kg], moisture [-], target_moisture [-], density [kg/m3]
Feeds_props = pd.read_excel(plib.Path(in_path, 'HTC_yield_HHV.xlsx'),
                            sheet_name='FeedsProperties', engine='openpyxl')
HC_yields_all = pd.read_excel(plib.Path(in_path, 'HTC_yield_HHV.xlsx'),
                              sheet_name='Yield_HC', engine='openpyxl') # [-]
HC_HHVs_all = pd.read_excel(plib.Path(in_path, 'HTC_yield_HHV.xlsx'),
                            sheet_name='HHV_HC', engine='openpyxl')  # [MJ/kg]
Gas_yields_all = pd.read_excel(plib.Path(in_path, 'HTC_yield_HHV.xlsx'),
                            sheet_name='Yield_Gas', engine='openpyxl')  # [-]
# Utility_factors_all: [kWh/kg feed]
Utility_factors_all = pd.read_excel(plib.Path(in_path, 'HTC_yield_HHV.xlsx'),
                                    sheet_name='UtilityFactors', engine='openpyxl') 
# PlantCapacity [t (dry)/y], Capex [$] 
HTCplantCapex = pd.read_excel(plib.Path(in_path, 'HTCplantCapex.xlsx'), 
                   sheet_name='WetCapacity', engine='openpyxl')

CapexFunc = CapexCurve(HTCplantCapex, fit_type='modified_econ_scale')  # [$]

n = 100  # grid resolution (nxn).  # use >= 30 for debug and >= 100 for results
# longitude, vertical lines, so horizontal scale
longs = np.linspace(-80, -72, n).repeat(n)  # creating all couples of coords
# latitude, horiontal lines so the vertical scale
lats = np.tile(np.linspace(40, 45, n), n)  # creating all couples of coords
df_latslongs = pd.DataFrame([lats, longs]).T  # put coords in pd.DF
df_latslongs.columns = ['lat', 'long']
# create a gpd.GDF with locations and geometries
geolatslongs = gpd.GeoDataFrame(df_latslongs,
                                geometry=gpd.points_from_xy(df_latslongs.long,
                                                            df_latslongs.lat))
# retrieve NY boundaries as a geometry
states = gpd.read_file(plib.Path(in_path, 'data/usa-states-census-2014.shp'))
# same coordinate system for map and locations: Spherical Mercator- GMaps
geolatslongs.crs = "EPSG:3395"
states.crs = "EPSG:3395"
NY = states[states['NAME'] == 'New York']
# select only locations inside NY to only run computations on those
locsInNY = gpd.tools.sjoin(geolatslongs, NY, how='inner')
locsLeftNY = gpd.tools.sjoin(geolatslongs, NY, how='left')
latslongs = list(zip(locsInNY.lat, locsInNY.long))  # create list of tuples


if 1:# food waste map
    FoodWasteDistribution(Spatialdata, feed_names, Q_max=40, factor=250)
    
#%%
# variables that are investigated as the final results in the following analyses
vrbls = [
    vrbl('En_gen', 'En. Gen. [GJ/day]', [100, 700], conv_factor=1e-3,
         annotate_lttr='a'),
    vrbl('NPV', 'NPV [M\$]', [-10, 20], conv_factor=1e-6),   
    vrbl('Carbon_Price', 'Carbon Price [\$/t$_{CO_{2, eq}}$]', [0, 250],
         pm_sign=[True, False], annotate_lttr='d'),
    vrbl('MinSellP_HC_eff', 'MSP$_{HC, eff}$ [\$/t]', [-50, 150],
         annotate_lttr='c'),
    vrbl('Net_CO2_emiss', 'Net CO$_2$ Emiss. [kt$_{CO_{2, eq}}$/y]', 
         [-14, -6], conv_factor=1e-3, annotate_lttr='b'),
    vrbl('CI_HC_mass', 'CI$_{HC}$ [t$_{CO_{2, eq}}$/t$_{HC}$]'),
    # vrbl('CI_HC_energy', 'CI$_{HC}$ [g$_{CO_{2, eq}}$/MJ]',[])
    ]
if 1:  # base case for paper (conditions are given by the parameters)
    subfolder='Case_Q2_r40_T220_t1'       
    prmtrs = [prmtr('Q_min', 1.81437, 'Q$_{min}$', 't/week'),
              prmtr('buf_radius', 40.2235, 'r$_{buffer}$', 'km'),
              prmtr('T_HTC', 220, 'T$_{HTC}$', '°C'),
              prmtr('t_HTC', 1, 't$_{HTC}$', 'h')]
    # general heat maps for spatial distr of each vrbls using the prmts
    HM_2_40_220_1 = MultiLoc(latslongs, prmtrs=prmtrs, vrbls=vrbls,
                             subfolder=subfolder)
    # selection of n best locations that maxim (or minimize) choic_var 
    BL_2_40_220_1_NPV = BestLocSelect(HM_2_40_220_1, prmtrs, vrbls,
                                      choice_var='NPV', subfolder=subfolder)                          
    BL_2_40_220_1_En_gen = BestLocSelect(HM_2_40_220_1, prmtrs, vrbls,
                                          choice_var='En_gen',
                                          subfolder=subfolder)
    BL_2_40_220_1_Net_CO2_emiss = BestLocSelect(HM_2_40_220_1, prmtrs,
                                                vrbls, choice_var='Net_CO2_emiss',
                                                choice_maximize=False,
                                                subfolder=subfolder)  
#%%

if 1 :  # fitting for HTC plant Capex plot
    PlCap_conv = 1/365  
    Cap_conv = 1e-6
    x = HTCplantCapex['Capacity_t_y'].values  # t
    y = HTCplantCapex['Capex_USD'].values  # $
    xmax = int(np.max(np.concatenate([BL_2_40_220_1_Net_CO2_emiss.PlantCapacity.values,
                             BL_2_40_220_1_NPV.PlantCapacity.values, x])))*2
    labs = ['y = p$_1$ $\cdot$x + p$_2$', 'y = p$_1$$\cdot$(x/p$_2$)$^{0.6}$',
            'y = p$_1$ + (x/p$_2$)$^{0.6}$']
    fig, ax, axt, fig_par = FigCreate(1, 1, 0, paper_col=.78, hgt_mltp=1)
    ax[0].plot(x*PlCap_conv, y*Cap_conv, color='k', marker=mrkrs[0], linestyle='None',
               label='Literature data')
    ax[0].plot(BL_2_40_220_1_Net_CO2_emiss.PlantCapacity*PlCap_conv,
               BL_2_40_220_1_Net_CO2_emiss.Capex*Cap_conv, color=clrs[0],
               marker=mrkrs[1], linestyle='None', label='Opt. Loc. (CO$_2$)')
    ax[0].plot(BL_2_40_220_1_NPV.PlantCapacity*PlCap_conv,
               BL_2_40_220_1_NPV.Capex*Cap_conv, color=clrs[1],
               marker=mrkrs[2], linestyle='None', label='Opt. Loc. (NPV)') 
    print('Fitting parameters [p$_1$, p$_2$]: ')
    for g, fit in enumerate(['linear', 'econ_scale',
                             'modified_econ_scale']):
        fun, fun_prmtrs = CapexCurve(HTCplantCapex, fit_type=fit,
                                     print_prmtrs=True)
        
        print(fit, fun_prmtrs)
        ax[0].plot(np.linspace(0, xmax, 100)*PlCap_conv,
                   fun(np.linspace(0, xmax, 100))*Cap_conv, 
                   color=clrs[2+g], linestyle=lnstls[g], label=labs[g])
    FigSave('CapexFit', out_path, fig, ax, axt, fig_par,
            xLim=[0, xmax*PlCap_conv], 
            yLim=[0, 20],
            xLab='Plant Capacity [t$_{wet}$/day]', legend='best', ncol_leg=1,
            yLab='Capex [M$]', tight_layout=True)  
#%%  CumulSystemPerf
if 1: # analysis of the cumulative effect of multiple plants, using the first
    # n_best locations for NYS. Base case prmtrs are used
    prmtrs = [prmtr('Q_min', 1.81437, 'Q$_{min}$', 't/week'),
              prmtr('buf_radius', 40.2235, 'r$_{buffer}$', 'km'),
              prmtr('T_HTC', 220, 'T$_{HTC}$', '°C'),
              prmtr('t_HTC', 1, 't$_{HTC}$', 'h')]

    vrbls_CSP = [vrbl('NPV', 'NPV [M\$]',[-15, 20], conv_factor=1e-6,
                      label_no_unit='NPV'), 
                 # the conversion from emiss to reduction is in the conv_factor
                 vrbl('Net_CO2_emiss', 'CO$_2$ Em. Reduc. [kt$_{CO_{2, eq}}$/y]', #[-20, -4],
                      conv_factor=-1*1e-3, label_no_unit='CO$_2$ Em. Reduc.')]
    CumulSystemPerf(HM_2_40_220_1, prmtrs, vrbls_CSP,
                    choice_var='Net_CO2_emiss',
                    choice_maximize=False, n_best=13,
                    yLim=[-160, 0], ytLim=[0, 60],
                    xTicks=[1, 5, 9, 13], 
                    yTicks=[-150, -100, -50, 0],
                    ytTicks=[0, 20, 40, 60])
    vrbls_CSP = [vrbl('NPV', 'NPV [M\$]',[-15, 20], conv_factor=1e-6,
                      label_no_unit='NPV'), 
                 # the conversion from emiss to reduction is in the conv_factor
                 vrbl('En_gen', 'En. Gen. [TJ/day]', [0, 700], conv_factor=1e-6,
                      label_no_unit='En. Gen.')]
    CumulSystemPerf(HM_2_40_220_1, prmtrs, vrbls_CSP,
                    choice_var='Net_CO2_emiss', choice_maximize=False,  n_best=13,
                    yLim=[-160, 0], ytLim=[0, 4],
                    xTicks=[1, 5, 9, 13], 
                    yTicks=[-150, -100, -50, 0],
                    ytTicks=[0, 1, 2, 3, 4])
#%%    
if 0: # example of an additional case at different conditions (not in the paper)
    subfolder='Case_Q2_r40_T250_t1'       
    prmtrs = [prmtr('Q_min', 1.81437, 'Q$_{min}$', 't/week'),
              prmtr('buf_radius', 40.2235, 'r$_{buffer}$', 'km'),
              prmtr('T_HTC', 250, 'T$_{HTC}$', '°C'),
              prmtr('t_HTC', 1, 't$_{HTC}$', 'h')]
    HM_2_40_250_1 = MultiLoc(latslongs, prmtrs=prmtrs, vrbls=vrbls,
                             subfolder=subfolder)   
    BL_2_40_250_1_NPV = BestLocSelect(HM_2_40_250_1, prmtrs, vrbls,
                                      choice_var='NPV', subfolder=subfolder)                          
    BL_2_40_250_1_En_gen = BestLocSelect(HM_2_40_250_1, prmtrs, vrbls,
                                         choice_var='En_gen', subfolder=subfolder)
    BL_2_40_250_1_Net_CO2_emiss = BestLocSelect(HM_2_40_250_1, prmtrs, vrbls,
                                                choice_var='Net_CO2_emiss',
                                                choice_maximize=False,
                                                subfolder=subfolder)
#%%
if 1:
    # single location base case cost and emission break down

    BL_labels = ['NYC', 'Buffalo', 'Rochester','Syracuse', 'Albany', 'LastOne'] 
 
    for p in range(BL_2_40_220_1_NPV.index[-1] + 1):  # all but NYC
        if p == 0:  # if NYC different limits are needed
            xLim_cost_emiss=[[-50, 50], [-80, 80]]
            xTicks_cost_emiss=[[-50, -25, 0, 25, 50],
                               [-80, -40, 0 , 40, 80]]
        else:
            xLim_cost_emiss=[[-16, 16], [-24, 24]]
            xTicks_cost_emiss=[[-16, -8, 0, 8, 16],
                               [-24, -12, 0 , 12, 24]] 
        SingleLoc(feed_names=feed_names, Spatialdata=Spatialdata,
                  Feeds_props=Feeds_props, HC_yields_all=HC_yields_all,
                  HC_HHVs_all=HC_HHVs_all, Gas_yields_all=Gas_yields_all,
                  Utility_factors_all=Utility_factors_all, CapexFunc=CapexFunc,
                  latlong=(BL_2_40_220_1_NPV.lat[p], BL_2_40_220_1_NPV.long[p]),
                  cost_emiss_brkdwn=True, file_name=BL_labels[p],
                  xLim_cost_emiss=xLim_cost_emiss, 
                  xTicks_cost_emiss=xTicks_cost_emiss, report=True)   
#%%   
if 1: # sensitivity analysis
    # HTC conditions are also varied,
    sens_prmtrs2 = [
        sens_prmt(name='latlong', label='location*', unit='deg',
                  base_value=(43.08080808080808, -77.57575757575758), # Rochester
                  do_not_vary=True, n_steps=5),
        sens_prmt(name='T_t_HTC', label=r'HTC conditions**', unit='C',
                  steps=[(190, 1), (190, 3), (220, 1), (220, 3), (250, 1)]),
        sens_prmt(name='Q_min', label=r'Q$_{min}$', unit='t/week',
                  base_value=1.81437, positive=False),
        sens_prmt(name='buf_radius', label=r'r$_{buffer}$', unit='km',
                  base_value=40.2335),
        sens_prmt(name='transp_price_h', label=r'price$_{transp.\ liq}$',
                  unit='check', base_value=100, positive=False),
        sens_prmt(name='transp_price_d', label=r'price$_{transp.\ sol}$',
                  unit='check', base_value=2.50, positive=False),
        sens_prmt(name='price_HC_USD_t', label=r'HC$_{sell.\ price}$',
                  unit='$/t', base_value=16.5),
        sens_prmt(name='price_elect', label=r'price$_{electr.}$',
                  unit='$/kWh', base_value=0.055, positive=False),
        sens_prmt(name='price_NatGas', label=r'price$_{nat.\ gas}$',
                  unit='$/ft3', base_value=0.00673, positive=False),
        sens_prmt(name='labor_rate', label=r'labor rate',
                  unit='$/hr', base_value=27.5, positive=False),
        sens_prmt(name='P_water_treatment', label=r'price$_{water\ treat.}$',
                  unit='$/m$^3$', base_value=0.86, positive=False),
        sens_prmt(name='tipping_fee', label=r'tipping fee',
                  unit='$/t', base_value=78.264),
        sens_prmt(name='dscnt_rt', label=r'r$_{discount}$',
                  unit='-', base_value=0.07, positive=False)]
    sens_vrbls2 = [vrbl('NPV', '$\Delta$ NPV [M\$]', conv_factor=1e-6, label_no_unit='NPV'),
                  vrbl('Net_CO2_emiss', 
                       '$\Delta$ Net CO$_2$ Emiss. [kt$_{CO_{2, eq}}$/y]',
                       conv_factor=1e-3, 
                       label_no_unit='Net CO$_2$ Emiss.', to_maximize=False)]
          
    SensAnalysis(sens_prmtrs2, sens_vrbls2, 
                 file_name='LinAndNonLinSA',
                 legend=['lower right','upper right'], ncol_leg=2,
                 yLim=[[-20, 20], [-10, 10]],
                 yTicks=[[-20, -10, 0, 10, 20], [-10, -5, 0, 5, 10]])
    


#%%
# ECONOMY of scale 
if 1: # sensitivity analysis where non-linear values (e.g., location, 
    # HTC conditions) are varied, results are not directly comparable but
    # provide an idea of the extent of the effect of these parameters
    sens_prmtrs2 = [
        sens_prmt(name='latlong', label='location*', unit='deg',
                  base_value=(42.82828282828283, -73.77777777777777), 
                  do_not_vary=True),
        sens_prmt(name='T_t_HTC', label=r'HTC conditions*', unit='C',
                  base_value=(220, 1), do_not_vary=True),
        sens_prmt(name='Q_min', label=r'Q$_{min}$', unit='t/week',
                  base_value=2, positive=False,
                  do_not_vary=True),
        sens_prmt(name='buf_radius', label=r'r$_{buffer}$', unit='km',
                  base_value=40, steps=np.arange(0, 600, 5)),
        sens_prmt(name='dscnt_rt', label=r'r$_{discount}$',
                  unit='-', base_value=0.07, positive=False,
                  do_not_vary=True),
        sens_prmt(name='tipping_fee', label=r'tipping fee',
                  unit='$/t', base_value=71,
                  do_not_vary=True),
        sens_prmt(name='transp_price_d', label=r'price$_{transp.\ sol}$',
                  unit='check', base_value=2.50, positive=False,
                  do_not_vary=True),
        sens_prmt(name='transp_price_h', label=r'price$_{transp.\ liq}$',
                  unit='check', base_value=100, positive=False,
                  do_not_vary=True), 
        sens_prmt(name='price_HC_USD_t', label=r'HC$_{sell.\ price}$',
                  unit='$/t', base_value=16.5, do_not_vary=True),
        sens_prmt(name='price_elect', label=r'price$_{electr.}$',
                  unit='$/kWh', base_value=0.055, positive=False,
                  do_not_vary=True),
        sens_prmt(name='price_NatGas', label=r'price$_{nat.\ gas}$',
                  unit='$/kWh', base_value=0.00673, positive=False,
                  do_not_vary=True),
        sens_prmt(name='labor_rate', label=r'labor rate',
                  unit='$/hr', base_value=27.5, positive=False,
                  do_not_vary=True),
        sens_prmt(name='P_water_treatment', label=r'price$_{water\ treat.}$',
                  unit='$/m$^3$', base_value=0.86, positive=False,
                  do_not_vary=True)]    
    sens_vrbls2 = [vrbl('NPV', 'NPV [M\$]', conv_factor=1e-6, label_no_unit='NPV'),
                  vrbl('Net_CO2_emiss', 
                       'Net CO$_2$ Reduction [kt$_{CO_{2, eq}}$/y]',
                       conv_factor=-1*1e-3, 
                       label_no_unit='Net CO$_2$ Red.', to_maximize=False)]
    
    BL_labels = ['NYC', 'Buffalo', 'Rochester','Syracuse', 'Albany', 'LastOne']
    dfs_res = []
    for p in range(BL_2_40_220_1_NPV.index[-1] + 1): 
        sens_prmtrs2[0] = sens_prmt(name='latlong', label='location*', unit='deg',
                                   base_value=(BL_2_40_220_1_NPV.lat[p],
                                               BL_2_40_220_1_NPV.long[p]), 
                                   do_not_vary=True)     
        res, res_var = SensAnalysis(sens_prmtrs2, sens_vrbls2,
                                    file_name='EconScaleR_' + BL_labels[p],
                                    legend='center right', ncol_leg=1,
                                    yLim=[-20, 25], ytLim=[-100, 0],
                                    type_analysis='single_param',
                                    plt_prmt='buf_radius')
        dfs_res.append(res)
        
    file_name = 'EconomyOfScale_rbuff'     
    fig, ax, axt, fig_par = FigCreate(rows=2, cols=1, plot_type=0,
                                      paper_col=.78, hgt_mltp=1.5)
    BL_labels = ['Buffalo', 'Rochester','Syracuse', 'Albany', 'White Plains']
    for d, df in enumerate(dfs_res):
          col =  'buf_radius' 
          ax[0].plot(df[0].index, df[0][col], color=clrs[d], 
               linestyle=lnstls[d], label=BL_labels[d])
          ax[1].plot(df[1].index, df[1][col], color=clrs[d], 
                    linestyle=lnstls[d], label=BL_labels[d])            
    FigSave(file_name, out_path, fig, ax, axt, fig_par,
            xLim=[0, 600], xLab='r$_{buffer}$ [km]',
            yLim=[[-20, 20], [0, 100]], 
            yLab=['NPV [M$]', 'Net CO$_2$ Reduction [kt$_{CO_{2, eq}}$/y]'], 
            yTicks=[[-20, -10, 0, 10, 20], [0, 20, 40, 60, 80, 100]],
                # ytTicks=yTicks, 
                # subfolder=subfolder,
                tight_layout=True,
                ncol_leg=1, 
                legend='best')    
#%%
# ECONOMY of scale 
if 1: # sensitivity analysis where non-linear values (e.g., location, 
    # HTC conditions) are varied, results are not directly comparable but
    # provide an idea of the extent of the effect of these parameters
    sens_prmtrs2 = [
        sens_prmt(name='latlong', label='location*', unit='deg',
                  base_value=(42.82828282828283, -73.77777777777777), 
                  do_not_vary=True),
        sens_prmt(name='T_t_HTC', label=r'HTC conditions*', unit='C',
                  base_value=(220, 1), do_not_vary=True),
        sens_prmt(name='Q_min', label=r'Q$_{min}$', unit='t/week',
                  base_value=2, positive=False,
                  steps=np.arange(0, 30, .1)),
        sens_prmt(name='buf_radius', label=r'r$_{buffer}$', unit='km',
                  base_value=40,do_not_vary=True) ,
        sens_prmt(name='dscnt_rt', label=r'r$_{discount}$',
                  unit='-', base_value=0.07, positive=False,
                  do_not_vary=True),
        sens_prmt(name='tipping_fee', label=r'tipping fee',
                  unit='$/t', base_value=71,
                  do_not_vary=True),
        sens_prmt(name='transp_price_d', label=r'price$_{transp.\ sol}$',
                  unit='check', base_value=2.50, positive=False,
                  do_not_vary=True),
        sens_prmt(name='transp_price_h', label=r'price$_{transp.\ liq}$',
                  unit='check', base_value=100, positive=False,
                  do_not_vary=True), 
        sens_prmt(name='price_HC_USD_t', label=r'HC$_{sell.\ price}$',
                  unit='$/t', base_value=16.5, do_not_vary=True),
        sens_prmt(name='price_elect', label=r'price$_{electr.}$',
                  unit='$/kWh', base_value=0.055, positive=False,
                  do_not_vary=True),
        sens_prmt(name='price_NatGas', label=r'price$_{nat.\ gas}$',
                  unit='$/kWh', base_value=0.00673, positive=False,
                  do_not_vary=True),
        sens_prmt(name='labor_rate', label=r'labor rate',
                  unit='$/hr', base_value=27.5, positive=False,
                  do_not_vary=True),
        sens_prmt(name='P_water_treatment', label=r'price$_{water\ treat.}$',
                  unit='$/m$^3$', base_value=0.86, positive=False,
                  do_not_vary=True)]    
    sens_vrbls2 = [vrbl('NPV', 'NPV [M\$]', conv_factor=1e-6, label_no_unit='NPV'),
                  vrbl('Net_CO2_emiss', 
                       'Net CO$_2$ Reduction [kt$_{CO_{2, eq}}$/y]',
                       conv_factor=-1*1e-3, 
                       label_no_unit='Net CO$_2$ Red.', to_maximize=False)]
    
    BL_labels = ['NYC', 'Buffalo', 'Rochester','Syracuse', 'Albany', 'LastOne']
    dfs_res = []
    for p in range(BL_2_40_220_1_NPV.index[-1] + 1): 
        sens_prmtrs2[0] = sens_prmt(name='latlong', label='location*', unit='deg',
                                   base_value=(BL_2_40_220_1_NPV.lat[p],
                                               BL_2_40_220_1_NPV.long[p]), 
                                   do_not_vary=True)     
        res, res_var = SensAnalysis(sens_prmtrs2, sens_vrbls2,
                                    file_name='EconScaleR_' + BL_labels[p],
                                    legend='center right', ncol_leg=1,
                                    yLim=[-20, 25], ytLim=[-100, 0],
                                    type_analysis='single_param',
                                    plt_prmt='buf_radius')
        dfs_res.append(res)
    
    file_name = 'EconomyOfScale_Qmin'     
    fig, ax, axt, fig_par = FigCreate(rows=2, cols=1, plot_type=0,
                                      paper_col=.78, hgt_mltp=1.5)
    BL_labels = ['Buffalo', 'Rochester','Syracuse', 'Albany', 'White Plains']
    for d, df in enumerate(dfs_res):
          col =  'Q_min' 
          ax[0].plot(df[0].index, df[0][col], color=clrs[d], 
               linestyle=lnstls[d], label=BL_labels[d])
          ax[1].plot(df[1].index, df[1][col], color=clrs[d], 
                    linestyle=lnstls[d], label=BL_labels[d])            
    FigSave(file_name, out_path, fig, ax, axt, fig_par,
            xLim=[0, 30], xLab='Q$_{min}$ [t/week]',
            yLim=[[-15, 5], [0, 20]], 
            yLab=['NPV [M$]', 'Net CO$_2$ Reduction [kt$_{CO_{2, eq}}$/y]'], 
            yTicks=[[-15, -10, -5, 0, 5], [0, 5, 10, 15, 20]],
                # ytTicks=yTicks, 
                # subfolder=subfolder,
                tight_layout=True,
                ncol_leg=1, 
                legend='best') 

# %% Grafical Abstract and Cover Art graphics
if 1:
    subfolder = 'GraphAbst_CovArt'
    plib.Path(plib.Path(out_path, subfolder)).mkdir(parents=True, exist_ok=True)
    
    # plot for cover art using best case best locations
    BestLocSelect(HM_2_40_220_1, prmtrs=None, name='CoverArt_BL_CO2',
                  vrbls=vrbls, choice_var='Net_CO2_emiss',
                  choice_maximize=False, subfolder=subfolder,
                  transparency=True, negative=True,
                  cities = ['White\nPlains',
                            'Buffalo', 'Syracuse', 'Rochester', 'Albany'], 
                  lat_cities = [41.9, 42.5, 43.4, 42.8, 43.4], 
                  long_cities = [-74.7, -78.7, -75.7, -77.6, -74.3],
                  legend=False, df_best_to_csv=False) 
#%%
    BestLocSelect(HM_2_40_220_1, prmtrs=None, name='CoverArt_BL_CO2',
                  vrbls=vrbls, choice_var='Net_CO2_emiss',
                  choice_maximize=False, subfolder=subfolder,
                  transparency=False, negative=False,
                  cities = ['White\nPlains',
                            'Buffalo', 'Syracuse', 'Rochester', 'Albany'], 
                  lat_cities = [41.9, 42.5, 43.4, 42.8, 43.4], 
                  long_cities = [-74.7, -78.7, -75.7, -77.6, -74.3],
                  legend=True, df_best_to_csv=False) 
#%% Clean Map of NYS, all black 
    fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=0)
    ax[0] = NY.plot(color='k', figsize=(8, 8), aspect='1.4',
                    facecolor='k', edgecolor='k')
    ax[0].set_axis_off()   
    FigSave('NYS_k', out_path, fig, ax, axt, fig_par,
            subfolder=subfolder, tight_layout=True, transparency=True) 
    # %%  Clean Map of NYS, only border
    fig, ax, axt, fig_par = FigCreate(rows=1, cols=1, plot_type=0)
    ax[0] = NY.boundary.plot(color='k', figsize=(8, 8), aspect='1.4')
    ax[0].set_axis_off()
    FigSave('NYS', out_path, fig, ax, axt, fig_par,
            subfolder=subfolder, tight_layout=True, transparency=True) 

#%%   
    
    #%%
    vrbls_CSP = [vrbl('NPV', 'NPV [M\$]',[-20, 0], conv_factor=1e-6,
                      label_no_unit='NPV'), 
                 # the conversion from emiss to reduction is in the conv_factor
                 vrbl('Net_CO2_emiss', 'CO$_2$ Em. Reduc. [kt$_{CO_{2, eq}}$/y]', 
                      [-20, -4], conv_factor=-1*1e-3,
                      label_no_unit='CO$_2$ Em. Reduc.')]
    CumulSystemPerf(HM_2_40_220_1, prmtrs, vrbls_CSP,
                    choice_var='NPV', n_best=5,
                    subfolder='GraphAbst_CovArt', yLim=[-40, 0], ytLim=[0, 40],
                    xTicks=[1, 3, 5], yTicks=[-40, -20, 0], ytTicks=[0, 20, 40], 
                    axis_labels=False, paper_col=0.35, legend=False)
    
    #%% 
    
       

    


#%%
# =============================================================================
#  #%%
# if 0:  # sensitivity analysis 
#     # parameters to be varied for the sensitivity analysis    
#     sens_prmtrs = [
#         sens_prmt(name='latlong', label='location*', unit='deg',
#                   base_value=(43.08080808080808, -77.57575757575758), # Rochester
#                   do_not_vary=True, n_steps=5),
#         sens_prmt(name='T_t_HTC', label=r'HTC conditions*', unit='C',
#                   base_value=(220, 1), do_not_vary=True, n_steps=5),
#         sens_prmt(name='Q_min', label=r'Q$_{min}$', unit='t/week',
#                   base_value=1.81437, positive=False),
#         sens_prmt(name='buf_radius', label=r'r$_{buffer}$', unit='km',
#                   base_value=40.2335),
#         sens_prmt(name='transp_price_h', label=r'price$_{transp.\ liq}$',
#                   unit='check', base_value=100, positive=False),
#         sens_prmt(name='transp_price_d', label=r'price$_{transp.\ sol}$',
#                   unit='check', base_value=2.50, positive=False),
#         sens_prmt(name='price_HC_USD_t', label=r'HC$_{sell.\ price}$',
#                   unit='$/t', base_value=16.5),
#         sens_prmt(name='price_elect', label=r'price$_{electr.}$',
#                   unit='$/kWh', base_value=0.055, positive=False),
#         sens_prmt(name='price_NatGas', label=r'price$_{nat.\ gas}$',
#                   unit='$/ft3', base_value=0.00673, positive=False),
#         sens_prmt(name='labor_rate', label=r'labor rate',
#                   unit='$/hr', base_value=27.5, positive=False),
#         sens_prmt(name='P_water_treatment', label=r'price$_{water\ treat.}$',
#                   unit='$/m$^3$', base_value=0.86, positive=False),
#         sens_prmt(name='tipping_fee', label=r'tipping fee',
#                   unit='$/t', base_value=78.264),
#         sens_prmt(name='dscnt_rt', label=r'r$_{discount}$',
#                   unit='-', base_value=0.07, positive=False)]
#     # target variables for the sensitivity analysis
#     sens_vrbls = [vrbl('NPV', 'NPV [M$]', conv_factor=1e-6, label_no_unit='NPV'),
#                   vrbl('Net_CO2_emiss', 'Net CO$_2$ Emiss. [kt$_{CO_{2, eq}}$/y]', 
#                        label_no_unit='Net CO$_2$ Emiss.',
#                        conv_factor=1e-3, to_maximize=False)]
#     # sensitivity analysis for Rochester (prmts varied by +-0.5 adn 0.8)
#     BL_labels = ['NYC', 'Buffalo', 'Rochester','Syracuse', 'Albany', 'LastOne'] 
#     for p in range(BL_2_40_220_1_NPV.index[-1] + 1): 
#         sens_prmtrs[0] = sens_prmt(name='latlong', label='location*', unit='deg',
#                                    base_value=(BL_2_40_220_1_NPV.lat[p],
#                                                BL_2_40_220_1_NPV.long[p]), 
#                                    do_not_vary=True, n_steps=5)
#         SensAnalysis(sens_prmtrs, sens_vrbls, 
#                      file_name='SA220C1h_' + BL_labels[p],
#                      yLim=[[-20, 20], [-8, 8]],
#                      yTicks=[[-20, -10, 0, 10, 20], [-8, -4, 0, 4, 8]])
# =============================================================================